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I.  INTRODUCTION 


One  concept  which  could  increase  the  firepower  and  survivability  of 
the  "fire-and-forget"  class  of  helicopter  borne  missile  systems  would  be 
to  acquire  and  hand  off  simultaneously,  or  almost  so,  multiple  targets 
from  a  precision  pointing  and  tracxing  system  (PTS)  to  several  missile 
seekers  and  fire  in  a  very  short  period  of  time.  A  typical  over-all  fire 
control  configuration  is  shown  in  Figure  1-1.  The  pointing  and  tracking 
system  typically  consists  of  an  optics  train,  line  of  sight  (l OS)  stabi¬ 
lization  system,  forward  looking  inf'ared  (FUR)  imaging  system,  manual 
and  autotrack  system,  laser  range  finder  and  associated  electronics.  An 
imaging  missile  seeker  could  be  an  infrared  type.  It  is  assumed  that 
during  preflight  checkout  or  during  the  actual  flight,  the  lines  of  sight 
of  all  the  missile  seekers  are  aligned  with  the  line  of  sight  of  the  PTS. 
However,  due  to  gyro  drift,  boresighting  inaccuracies,  vehicle  vibration 
and  flexure,  etc.,  the  seekers  will  not  remain  tcresighted  with  the  PTS. 
Since  the  PTS  has  a  larger  field  of  view  { FOV )  in  both  axes  than  the  mis¬ 
sile  seekers,  it  is  expected  that  the  FOV  of  all  the  missile  seekers  will 
be  located  within  the  FOV  of  the  FTS,  The  multiple  target  problem  then 
becomes  that  of  locating  missile  seeker  LOS  aim  points  within  the  PTS 
field  of  view,  deciding  which  target  is  to  be  assigned  to  each  missile, 
generating  error  signals  to  the  torquers  in  order  to  slew  each  missile 
LOS  such  that  its  assigned  target  is  in  the  center  of  its  FOV,  and  initi¬ 
ating  automatic  seeker  tracking.  The  work  reported  in  this  final  report 
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assumes  that  some  method  is  available  for  acquiring,  recognizing,  and 
marking  potential  targets  within  the  PL  I R  FOV  (eg.,  MTI  radar  and/or  a 
target  recognizer) , 


Figure  1-1.  Helicopter  missile  fire  control  system. 

In  general,  the  task  of  locating  a  given  smaller  image  within  a 
larger  image  is  known  as  "image  registration".  The  smaller  image  is  re¬ 
ferred  to  as  the  window  or  the  reference  and  the  larger  image  is  called 
the  search  area.  Using  the  above  notation  with  the  multiple  target  prob¬ 
lem,  the  low  resolution  image  obtained  from  the  PTS  sensor  is  the  search 
area  and  the  high  resolution  image  obtained  from  each  missile  seeks',  is  a 
window.  Therefore,  ft.ere  is  one  search  area  and  more  than  one,  say  n, 
windows.  It  is  ass'. med  that  all  the  n  windows  are  completely  located 
within  the  search  area.  How  the  problem  of  multiple  image  registration 
can  be  defined  as  that  of  finding  n  subimages  within  the  search  area 
which  best  match  the  n  windows.  The  n  windows  W^,  W^,  . „ .  Wp  and  the 
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search  area  S  are  shown  in  Figure  1-2.  S  is  m  MxN  array  of  digital  pic¬ 
ture  elements  which  may  assume  one  of  G  possible  levels  on  the  gray  scale. 

0  i  S(i  ,.j )  <  G-l  ,  (1-1) 

for  1  <_  i  <_  M,  1  <_  j  5  N  . 

is  a  KxL  array  of  pixels  having  the  same  gray  scale  range. 

0  <_  W^U.m)  <  G-l  ,  (1-2) 

for  1  <_  a  <_  K,  1  <  m  <  L  and  k  =  1 ,  2,  ....  n  . 

Each  KxL  subimage  of  S  can  be  uniquely  identified  by  its  upper  left  cor¬ 
ner's  coordinates.  Let  S.  .  denote  the  KxL  subimage  of  S  whose  upper 
left  corner  is  (i ,  j ) . 

S.  ,U,m)  =  S(i+t-l  ,  j+m-1  )  ,  (1-3) 

*  »J 


for 

1  <  &  <  K,  1  <  ni 

<  L 

and 

1  <  i  <  M-K+l ,  1 

£  j 

If  S  and  W.  do  not  differ  in  pixel  resolution  and  rotation,  the  multiple 
image  registration  problem  reduces  to  that  of  finding  (i*,  j*)  such  that 
S.*  .*  best  matches  W,  ,  for  k  =  1,  2,  ...,  n. 

VJk  K 


Figure  1-2.  Search  area  and  windows. 
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Several  template  matching  and  feature  matching  methods  of  accomplish¬ 
ing  digital  image  registration  are  presented  in  detail  in  the  Final  Tech¬ 
nical  Report  of  Phase  I  of  contract  DAAH01  -80-00258.  It  has  been  shown 
that,  in  general,  feature  matching  algorithms  require  fewer  arithmetic 
operations  for  software  implementation  and  less  hardware  for  real-time 
hardware  implementation  than  the  template  matching  methods  such  as  "cor¬ 
relation"  and  "sequential  similarity  detection  algorithms",  if  the  number 
of  windows  is  sufficiently  large.  In  Chapter  II,  the  results  of  simula¬ 
tion  of  the  important  registration  methods  in  [1]  are  given  and  their 
registration  accuracies  are  compared.  Two  new  feature  matching  methods, 
one  based  on  the  Haar  transform  and  the  other  based  on  adjacent  pixel 
differences  along  the  rows  and  columns  of  a  digital  image,  are  also  pre¬ 
sented  in  Chapter  II.  Possible  hardware  implementations  of  the  algori¬ 
thms  discussed  in  Chapter  II  are  given  in  Chapter  III. 

The  main  thrust  of  this  effort  has  been  to  determine  the  reliability 
and  accuracy  of  existing  feature  matching  methods  and  to  define  new  meth¬ 
ods.  It  was  felt  that  in  order  to  reduce  computation  time  and/or  hardware 
requirements,  the  actual  matching  process  should  be  based  on  a  few  domi¬ 
nant  features  rather  than  image  gray  level  or  edge  content.  The  hardware 
vs.  time  tradeoffs  for  several  sequential  correlation  image  matching  tech¬ 
niques  and  for  the  feature  matching  technique  are  discussed  in  Chapter  IV. 
The  final  technique  in  Chapter  IV  suggests  a  technology  program  for  devel¬ 
oping  combined  target  recognizer  and  feature  matching  handoff  hardware. 

The  conclusions  and  recommendations  are  given  in  Chapter  V. 
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II.  IMAGE  REGISTRATION  METHODS  AND  SIMULATION  RESULTS 


The  problem  of  registering  a  given  image  within  a  larger  image  uses 
techniques  fundamental  to  the  disciplines  of  image  processing  and  pattern 
recognition.  The  image  registration  algorithms  are  broadly  classified  as 
template  matching  algorithms  and  feature  matching  algorithms.  Template 
matching  methods  such  as  "correlation"  and  "sequential  similarity  detec¬ 
tion  algorithms"  are  widely  used  for  the  determination  of  local  similar¬ 
ity  between  two  images  in  the  single  image  registration  problem.  For' 
multiple  image  registration,  it  has  been  shown  that  feature  matching 
methods  are  computationally  more  efficient  than  template  matching  methods 
if  the  number  of  windows  is  sufficiently  large  [1].  This  is  because  the 
amount  of  computation  for  template  matching  algorithms  is  directly  pro¬ 
portional  to  the  number  of  windows  whereas  in  feature  matching  algorithms 
it  is  not.  In  feature  matcl  ng  algorithms, features  are  extracted  for  all 
subimages  of  the  search  a ~ea  and  windows  only  once  and  the  matching  pro¬ 
cedure  is  repeated  once  for  each  window.  Since  the  amount  of  computation 
required  to  match  the  features  is  negligible  compared  to  that  required  to 
compute  the  features,  the  computational  efficiency  of  the  feature  match¬ 
ing  methods  increases  with  the  number  of  windows. 

Three  possible  feature  matching  methods  of  accomplishing  multiple 
image  registration  based  on  image  moments,  intraset  and  interset  distances, 
and  correlation  of  adjacent  pixels  are  described  in  detail  in  the  Final 
Technical  Report  of  Phase  I  of  contract  DAAH01 -80-C-0258.  All  three 
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methods  were  shown  to  be  computationally  more  efficient  than  the  most 
widely  used  cross  correlation  method  if  the  number  of  windows  is  suffi¬ 
ciently  large.  However,  the  most  important  performance  measure  of  an  al¬ 
gorithm  is  its  registration  accuracy.  During  this  phase  of  the  contract, 
an  effort  was  made  to  determine  the  accuracy  of  the  above  feature  match¬ 
ing  methods  through  simulation.  In  addition  to  the  above  methods,  two 
new  feature  matching  methods,  one  based  on  the  Haar  transform  and  the 
other  based  on  adjacent  pixel  differences  along  the  rows  and  columns  of  : 
digital  image  are  also  presented  in  this  chapter.  The  performance  (reg¬ 
istration  reliability)  of  each  method  is  compared  with  that  of  a  binary 
correlation  algorithm. 

A.  Scenes  Used  for  Simulation 

Four  different  sets  of  images  for  targets  at  four  different  ranges 
are  used  for  simulation.  Each  set  consists  of  a  low  resolution  FLIR  image 
and  four  high  resolution  FLIR  images  (i.e.,  four  targets  at  each  range). 
The  fields  of  view  of  the  high  resolution  images  are  completely  located 
within  the  field  of  view  of  the  low  resolution  image.  The  ranges  to  the 
targets  located  at  the  four  different  positions  are  given  in  Table  1. 

Table  1.  Distance  between  the  targets  and  the  sensors. 


Position 

Di stance 

A 

3  km 

B 

3.5  km 

C 

4  km 

D 

5  km 

i 
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All  image  registration  methods  presented  in  this  report  are  based  on 
the  assumotion  that  the  search  area  and  windows  do  not  differ  in  rotation 
and  spatial  pixel  resolution.  However,  due  to  the  difference  in  s  .-nsor 
characteristics  the  low  resolution  FUR  image  (search  area)  and  high  res¬ 
olution  FLIR  images  (windows)  do  not  have  the  same  spatial  pixel  resolu¬ 
tion.  In  order  to  equalize  the  spatial  pixel  resolution,  the  images  are 
preprocessed  as  described  below. 

1.  Each  field  of  the  unprocessed  high  resolution  FLIR  image  is 
a  406  x  392  array  of  digital  picture  elements  which  can  as¬ 
sume  one  of  256  possible  levels  on  the  gray  scale.  The  hor¬ 
izontal  and  vertical  scale  factors  which  are  used  to  reduce 
the  pixel  resolution  of  high  resolution  FLIR  images  are  pre¬ 
determined  to  be  5.0  and  4.108..  respectively.  Reducing  the 
resolution  using  the  above  scale  factors  yields  an  array  of 
size  98  x  78. 

2.  Each  field  of  unprocessed  low  resolution  FLIR  imagery  con¬ 
tains  406  x  392  pixels  where  each  pixel  can  assume  one  of 
256  possible  gray  levels.  The  horizontal  and  vertical  scale 
factors  which  reduce  the  resolution  of  the  low  resolution 
FLIR  image  to  that  of  the  reduced  high  resolution  image  are 
predetermined  to  be  1.25  and  1.0,  respectively.  Reducing 
the  resolution  using  the  above  scale  factors  yields  a  re¬ 
duced  low  resolution  FLIR  image  of  size  406  x  313.  The  1.25 
horizontal  scale  is  necessary  in  order  to  reduce  the  effec¬ 
tive  sampling  rate  for  the  digitized  video  from  12.5  MHz  to 


lii  MHz. 
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The  above  preprocessed  images  are  inputs  to  each  of  the  image  registra¬ 
tion  algorithms  simulated. 

B.  Binary  Correlation  Algorithm 

Cross  correlation  is  the  most  widely  used  method  of  accomplishing 
digital  image  registration.  If  the  window  and  search  area  are  preproces¬ 
sed  and  transformed  into  binary  images  to  accomplish  image  registration, 
a  substantial  data  reduction  and  savings  in  computation  can  be  achieved. 
Use  of  binary  images  not  only  reduces  the  computational  burden  of  the 
correlation  algorithm  for  software  implementation,  but  also  reduces  the 
complexity  of  a  real-time  hardware  implementation.  A  correlation  algori¬ 
thm  which  is  easily  irnpleinentable  in  real-time  using  simple  digital  hard¬ 
ware  (exclusive-NOR  gates)  is  given  in  Reference  [2],  Simulation  results 
given  in  References  [3,  4]  give  the  accuracy  trade-offs  for  binary  corre¬ 
lation.  Because  of  the  previous  analysis  and  hardware  design,  the  binary 
correlator  is  used  as  a  basis  for  comparison  of  the  performance  of  fea¬ 
ture  matching  registration  algorithms  in  this  report.  The  registration 
process  using  binary  images  is  described  in  the  following  steps. 

Step  1:  A  KxL  reference  image  W  is  selected  from  the  reduced  high  reso¬ 
lution  FLIP,  image. 

Step  2:  The  mean  pixel  value  of  W  is  computed.  The  reference  is  then 
transformed  to  a  binary  image  by  quantizing  pixels  with  values 
greater  than  or  equal  to  to  ones  and  those  with  values  less 
than  to  zeros. 

Step  3:  Each  KxL  subimage  of  the  search  area  (preprocessed  low  resolution 
FLIR  image)  is  also  quantized  to  two  levels  in  the  same  manner 
as  the  reference  image. 
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Step  4:  The  binary  images  are  then  correlated.  The  correlation  function, 

R ( i , j ) ,  which  is  a  measure  of  similarity  between  W  and  S.  . ,  is 

i » J 

computed  as 

K  t 

R(i,j)  =  KL  -  E  z  IS.  . -  W(£,m))  (2-1) 

i=l  m=l  1,J 

for  1  ^  i  <_  M-K+l  and  1  <  j  £  N-L+l  . 

Maximum  similarity  is  expected  at  the  registration  point 
The  four  preprocessed  low  resolution  images  and  sixteen  high  resolu¬ 
tion  images  described  iri  a  previous  section  are  used  for  simulation.  For 
each  scene  used  for  simulation,  the  following  three  cases  are  considered. 
Case  1:  Reference  images  of  size  32  x  32  are  chosen  to  include  the  prom¬ 
inent  features  in  the  reduced  high  resolution  FLIR  images.  In 
other  words,  references  are  not  necessarily  taken  from  the  center 
of  the  field  of  view  of  the  reduced  high  resolution  image.  For 
this  case  the  range  values  given  in  Table  1  are  meaningless 
since  those  ranges  refer  to  the  target  in  the  center  of  the 
field  of  view. 

Case  ?.:  Reference  images  of  size  32  x  32  are  chosen  from  the  center  of 
the  fields  of  view  of  the  reduced  high  resolution  images. 

Case  3:  Reference  images  of  size  32  x  64  are  chosen  from  the  center  of 
the  fields  of  view  of  the  reduced  high  resolution  images.  For 
this  case,  one  of  the  reference  images  (position  A,  target  # 1) 
is  not  completely  located  within  the  search  area  and  therefore 
only  fifteen  scenes  are  simulated. 

Correlation  simulation  results  for  each  of  the  above  cases  are  given 
in  detail  iri  the  monthly  reports  of  December  1980  through  June  1981  and 
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are  not  given  in  this  report.  A  summary  of  results  is  presented  in  Table 
2.  In  Table  2,  a  "1"  indicates  that  the  true  registration  point  appeared 
as  the  first  maximum  in  the  correlation  surface.  Similarly,  a  "2",  "3", 
or  "4"  indicate  that  the  registration  point  appeared  as  second,  third  or 
fourth  maximum,  respectively.  An  "X"  indicates  that  the  registration 
point  was  not  within  the  first  four  maxima.  The  first  four  maxima  in  the 
correlation  surface  are  found  by  finding  the  maximum  value  first.  The 
second  maximum  is  found  by  masking  out  an  array  of  21  x  21  pixels  centered 
about  the  first  maximum  and  then  searching  the  remaining  correlation  sur¬ 
face  for  its  maximum  value.  Third  and  fourth  maxima  are  found  similarly. 
For  Case  1  the  true  registration  point  appeared  as  the  first  maximum  in 
9  out  of  16  cases,  was  within  the  first  four  maxima  for  3  additional  cases, 
and  was  not  within  the  first  four  maxima  for  the  remaining  4  cases.  These 
results  and  the  results  for  Cases  2  and  3  can  be  found  in  Table  3.  When 
reference  images  were  chosen  to  include  prominent  features  in  the  reduced 
high  resolution  image,  the  performance  of  the  binary  correlation  algorithm 
was  satisfactory.  For  the  other  two  cases,  the  method's  performance  was 
poor. 


C.  Image  Registration  Using  Correlation  of  Adjacent  Pixels 
A  method  of  accomplishing  digital  image  registration  based  on  corre¬ 
lation  of  adjacent  pixels  is  given  in  Reference  1.  In  Reference  1,  this 
method  was  shown  to  be  computationally  more  efficient  than  any  of  the 
template  matching  or  feature  matching  algorithms  studied  to-date.  Multi¬ 
plication  and  addition  requirements  for  this  method  and  other  important 
registration  methods  are  derived  in  [1],  The  number  of  additions 
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X 


X 
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Table  3.  Correlation  simulation  results. 

Registration  point  Registration  point  Registration  point 
is  first  maximum  2nd,  3rd,  or  4th  not  in  four  maxima 

maximum 

Case  19  3  4 

Case  23  1  12 

Case  32  2  11 

and  multiplications  required  by  various  registration  methods  as  a  func¬ 
tion  of  the  number  of  windows  of  size  32  x  32  within  a  search  area  of 
size  240  x  256  are  shown  in  Figures  2-1  and  2-2,  respectively.  The  equiv¬ 
alent  number  of  additions  required  by  each  method  as  a  function  of  the 
number  of  windows  is  given  in  Figure  2-3.  Figure  2-3  is  derived  from 
Figures  2-1  and  2-2  assuming  that  each  real  multiplication  is  equivalent 
to  three  real  additions.  From  the  above  figures  it  is  clear  that,  even 
for  single  image  registration,  the  method  based  on  correlation  of  adjacent 
pixels  requires  fewer  arithmetic  operations  than  the  standard  cross  cor¬ 
relation  algorithm.  Since  the  feature  vector  of  the  subimage  S.,,  .  or 

i+l  ,J 

S.  .  ,  can  be  computed  from  the  feature  vector  of  the  subimage  S.  .  with 

1  3 J+  '  1  » J 

very  few  arithmetic  operations,  this  method  is  very  promising  for  real¬ 
time  implementation.  In  order  to  determine  the  reliability  and  accuracy 
of  the  correlation  of  adjacent  pixels  method,  simulation  was  conducted 
both  with  and  without  intensity  level  normalization  of  images.  The  reg¬ 
istration  technique  and  simulation  results  are  presented  next  for  all 
cases  considered. 

Case  1:  In  Case  1,  an  attempt  is  made  to  accomplish  image  registration 
without  intensity  level  normalization  or  any  other  preprocessing 


5.0x10 


Figure  2-1.  Number  of  additions  to  register  n  windows 
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Figure  2-2.  Number  of  multiplications  to  register  n  windows. 


5.0x10 


Figure  2-3.  Number  of  equivalent  additions  to  register  n  windows. 
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Step  1 

Step  2 

Step  3 

Step  4 

Step  5 


of  i  ^es.  The  window  and  search  area  whose  elements  can  as¬ 
sume  one  of  256  levels  on  the  gray  scale  are  used  with  the  reg¬ 
istration  procedure  as  given  in  the  seven  steps  below. 

A  reference  image  W  of  size  KxL  is  taken  from  the  reduced  high 
resolution  FLIR  image  to  include  prominent  features  of  the  image. 
The  normalized  correlation  between  adjacent  pixels  of  each  row 
is  computed  as 

l-l 

l  W(«,m)  W(.e,m+1) 


X,  L  o 

i  Vr(s.,m) 
m-1 

for  i  -  1 ,  2,  3,  . . . ,  K  . 


The  normalized  correlation  between  the  adjacent  pixels  of  each 
column  is  computed  as, 


K-l 

Z  W  U  ,m )  W(t-il,m) 

a  If] _ 

°m  K  ? 

Z  W  ( «.  ,m) 

£=1 

formal,  2,  3,  ...,L  . 


(2-3) 


The  (K+L)-dimensional  feature  vector,  V^,  of  the  window  W  is 
given  by 

~  I- p  1 p 2  **'  pK°la2  **’  °L"^  (2-4) 


Feature  vectors  for  all  KxL  subimages  of  the  search  area  of  size 

MxN  are  computed.  Let  V.  .  be  the  feature  vector  of  the  subimage 

i ,3 

S.  .. 

i  ,J 
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Step  6:  Finally,  an  error  surface  of  size  (M-K+l )  x  (N-L-tl)  is  computed. 
The  elements  of  the  error  surface,  E,  are  defined  to  be 

E(i.j)  =  llVVi,jli2  ’  (2-5) 

for  1  <  i  _<  M-K+l  and  1  <  j  <  N~L+i  . 

E(i,j)  is  a  measure  of  dissimilarity  between  W  and  S.  .  where 

'  90 

minimum  dissimilarity  is  expected  at  the  registration  point. 

Step  7;  Four  minima  in  the  error  surface  are  found  by  finding  the  mini¬ 
mum  value  first.  The  second  minimum  is  found  by  masking  out  an 
array  of  21  x  21  pixels  centered  about  first  minimum  and  then 
searching  the  remaining  error  surface  for  its  minimum  value. 

Third  and  fourth  minima  are  found  similarly. 

Using  the  sixteen  pairs  of  reduced  low  and  high  resolution  FLIR  im¬ 
ages,  the  above  method  was  simulated  for  reference  images  of  size  32  x  32. 
The  method  failed  to  identify  the  registration  point  for  all  scenes  and 
therefore  simulation  results  for  this  case  are  not  included  in  this  re¬ 
port. 

Case  2:  In  order  to  obtain  meaningful  results  from  any  of  the  image  reg¬ 
istration  algorithms,  it  is  necessary  to  preprocess  the  window 
and  subimages  of  the  search  area  such  that  the  mean  and  standard 
deviation  of  their  pixel  values  are  equal.  This  is  called  inten¬ 
sity  level  normalization  and  is  normally  required  when  the  window 
and  the  search  area  are  obtained  from  different  sensors  with 
different  d.c.  gain  and  bias.  Ideally,  the  window  and  each  sub¬ 
image  of  the  search  area  are  preprocessed  to  have  zero  mean  and 
unity  standard  deviation.  In  this  section,  the  effect  of 
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intensity  level  normalization  of  the  window  and  subimages  of  the 
search  area  is  presented.  The  registration  steps  and  simulation 
results  are  given  below. 

Step  1:  A  reference  image  W  of  size  KxL  is  taken  from  the  reduced  high 
resolution  image  to  include  prominent  features  of  the  image. 

Step  2:  Let  be  the  mean  and  ow  be  the  standard  deviation  of  pixel 

values  of  W.  The  intensity  level  normalized  image  of  the  window, 
W ,  is  computed  as 

W(it,m)  -  u,. 


W‘  U,m)  - 


JW 


(2-6) 


JW 


for  1  <  t  <  K  and  1  <  m  <  L  . 


Step  3:  Since  the  window  is  normalized,  it  is  no  longer  necessary  to  com¬ 
pute  the  normalized  correlation  between  adjacent  pixels  of  rows 
and  columns  of  W:.  The  unnormalized  correlation  between  the  ad¬ 
jacent  pixels  of  each  row  is  computed  as 

L-l 


P  =  Z 
£  m=l 


W  U,m)  W'  U,m+1 ) 
for  i  =  1,  2,  3,  ...,  K  . 


Step  4:  The  unnormalized  correlation  between  the  adjacent  pixels  of  each 
column  is  computed  as 


=  l  W'  (t,m)  W'  ( 2,4 1  ,m) 

4=1 

for  m  =  1 ,  2,  3,  . . . ,  L  . 

The  remainder  of  the  image  registration  procedure  is  the  same  as  that  of 
Case  1.  Simulation  results  using  reference  images  of  size  32  x  32  are 
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presented  in  Table  4.  Frorn  Table  4  it  is  clear  that  the  true  registra¬ 
tion  point  appeared  as  the  first  minimum  for  9  out  of  16  scenes.  The 
registration  point  was  within  the  first  four  minima  for  four  additional 
scenes  and  was  not  within  first  four  minima  for  the  remaining  three  scenes 
The  performance  of  the  registration  method  based  ori  the  correlation  of 
adjacent  pixels  for  intensity  level  normalized  images  is  comparable  with 
that  of  the  binary  correlation  algorithm.  From  Case  1  and  Case  2  it  is 
concluded  that  intensity  level  normalization  of  images  improves  perfor¬ 
mance  of  the  method  and  that  it  is  not  possible  to  use  the  registration  al 
gorithm  based  on  correlation  of  adjacent  pixels  without  normalizing  the 
imayes.  However,  if  each  subimage  of  the  search  area  is  normalized  inde¬ 
pendently,  it  is  not  possible  to  compute  the  feature  vector  of  S.  .  .  or 
S.  ,  .  from  the  feature  vector  of  S.  .  as  given  in  Reference  1.  1  his  in- 

1  +  I  5  J  1  7  J 

creases  the  amount  of  computation  for  software  implementation  and  makes 
the  real-time  hardware  implementation  more  complex  and  less  attractive. 
Therefore,  an  effort  is  made  to  develop  an  image  preprocessing  algorithm 
which  has  almost  the  same  effect  on  registration  accuracy  as  that  of  in¬ 
tensity  level  normalization  of  images,  but  allows  computation  of  V.  ,  . 

i  >  »  j 

or  ^i  j+1  ^rom  ^i  j'  ®ne  suc^  PreProcess’'n9  algorithm  is  presented  next. 

Case  3:  A  possible  reason  for  the  unsatisfactory  performance  of  the 

image  registration  method  based  on  correlation  of  adjacent 

pixels  using  an  unnormalized  window  and  search  area  is  given 

below.  Let  W  and  S ..  ..  be  related  as 
i  »J 

W(£,m)  =  a  S.*  .*  U,m)  +  b  , 

1  » J 

for  1  <  fc  <  K  and  1  <  m  <  L  , 


(2-8) 
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Step  1 


where  a  and  b  are  constants.  The  above  is  due  to  the  difference 
in  gain  and  d.c.  bias  of  the  two  imaging  systems.  From  Equation 
(2-8)  it  is  clear  that  the  feature  vectors  of  W  and  S.*  .*  are 
not  equal  if  b  f  0..  However,  if  the  two  preprocessed  images  W' 
and  S'  ,*  are  related  by 

■  J  J 

W'U,m)  =  a  S'  (i,m)  ,  (2-9) 

*  » J 

for  1  <  t  <  K  and  1  £  m  £  L  , 


then  the  feature  vectors  of  W*  and  S'.*  .*  computed  as  described 
in  Case  1  will  be  equal.  A  method  of  accomplishing  the  above  is 
given  in  the  steps  below. 

The  value  of  each  pixel  of  the  reduced  high  resolution  FLIP 
image  is  replaced  by  its  value  minus  the  average  pixel  value  of 
the  K1  x  LI  array  centered  about  it.  This  process  can  be  better 
understood  from  Figure  2-4  where  X  is  the  reduced  high  resolu¬ 
tion  FLIR  image.  X ( i ,  j )  is  the  original  value  of  pixel  (i,j) 
and  the  new  value  X'(i,j)  is  computed  as 


X'(iJ)  =  X(i,j)  -  y.  . 

*  i  J 


where  y.  . 


1 


i,j  K1  x  1.1 


Kl-1 

2 


Ll_-X 

2 


(2-10) 


J]  X(i+k,  jn) 

LI  -1 


c  c 


and  where  K1  and  11  are  odd  numbers, 


Hew  values  for  all  other  pixels  are  computed  similarly.  If  the 
input  image  is  an  MxN  array,  the  preprocessed  image  will  be  an 
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(M-Kl+1)  x  (N-Ll+1)  array  of  valid  pixels.  Reference  image  W 
of  size  KxL  is  taken  from  the  above  (M-Kl+1)  x  (N-Ll+1)  array  of 
valid  pixels. 

Step  2:  The  reduced  low  resolution  FLIR  image  is  preprocessed  as  described 
in  Step  1  and  is  used  as  the  search  area  with  the  registration 
procedure  the  same  as  that  described  in  Case  1.  Simulation  re¬ 
sults  for  a  reference  image  of  size  32  x  32,  and  for  various 
values  of  K1  and  L.1  are  given  in  Table  4.  For  K1  =  LI  =11,  the 
true  registration  point  appeared  as  the  first  minimum  in  the  er¬ 
ror  surface  for  9  out  of  16  scenes  simulated.  Registration  point 
was  within  the  first  four  minima  for  four  additional  scenes  and 
was  not  within  the  first  four  minima  for  the  remaining  three 
scenes.  A  suumiary  of  trie  results  in  Table  4  can  be  found  in 
Table  6.  The  performance  of  the  method  for  Kl  =  LI  -  11  is  as 
good  as  that  of  binary  correlation.  Since  each  subimage  of  the 
search  area  is  not  preprocessed  separately  and  since  intensity 
level  normalization  is  no  longer  required  after  preprocessing, 

the  feature  vectors  V.,,  .  and  V.  ...  can  be  computed  from  V.  • 

i  +  *  >  J  i » J  "*'  *  ijJ 

with  very  few  arithmetic  operations  as  described  in  Reference  1. 

Case  4:  The  method  of  computing  feature  vectors  is  the  same  as  that  of 
Case  3  (Kl  =  LI  =  11).  However,  a  different  criterion  is  used 
for  matching  the  feature  vectors.  Instead  of  computing 
[  |V|,  -  V-  .  |  |  as  a  measure  of  dissimilarity  between  W  and  S..  . , 

w  1  }  j  i  >  j 


vT  V 

w  Y.i 


w 


|V.  . 


T(i,j)  = 


(2-11) 


23 


is  computed  as  a  measure  of  similarity  where  maximum  similarity 
is  expected  at  the  registration  point.  The  procedure  of  finding 
the  first  four  maxima  in  the  correlation  surface  is  the  same  as 
the  one  described  earlier.  The  simulation  results  for  the  six¬ 
teen  pairs  of  scenes  simulated  are  given  in  Table  4.  The  true 
registration  point  appeared  as  the  first  maximum  in  the  correla¬ 
tion  surface  for  12  out  of  16  scenes.  Registration  point  was 
within  the  first  four  maxima  for  two  additional  scenes  and  was 
not  within  the  first  four  maxima  for  the  two  remaining  scenes. 


Figure  2-4.  Layout  for  preprocessing  the  reduced 
high  resolution  image. 

Case  5:  In  Case  4,  each  reference  was  chosen  to  include  prominent  fea¬ 
tures  of  its  corresponding  preprocessed  reduced  high  resolution 
image.  In  otherwords,  the  reference  image  was  riot  necessarily 
taken  from  the  center  of  the  field  of  view  of  the  reduced  high 
resolution  image.  In  Case  5,  the  simulation  is  repeated  by 
choosing  the  reference  images  from  the  center  of  field  of  view. 
Simulation  results  for  reference  images  of  size  32  x  32  and 
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32  x  64  are  given  in  Table  5.  For  reference  of  size  32  x  32, 
the  registration  point  appeared  as  the  first  maximum  for  only  3 
out  of  16  scenes.  Registration  point  appeared  as  second,  third 
or  fourth  maximum  for  two  additional  scenes  and  was  not  within 
the  four  maxima  for  the  remaining  eleven  scenes.  Even  refer¬ 
ence  images  of  size  32  x  64  gave  similar  results  and  the  perfor¬ 
mance  was  not  improved. 

0.  Image  Registration  Based  on  Correlation  of 
Adjacent  Rows  and  Columns 

This  method  is  also  based  on  correlation  of  adjacent  pixels.  Instead 
of  computing  a  (K+L)-dimensional  feature  vector  for  each  KxL  array  based 
on  correlation  of  adjacent  pixels  in  each  row  and  column,  a  ( K+L-2 ) -di - 
mensional  feature  vector  is  computed  based  on  correlation  of  adjacent  rows 
and  adjacent  columns  as  described  in  the  following  steps. 

Step  1:  The  reduced  high  and  low  resolution  images  are  preprocessed  using 
an  algorithm  which  replaces  the  value  of  each  pixel  by  its  value 
minus  the  average  pixel  value  of  the  11  x  11  array  centered  about 
the  pixel . 

Step  2:  The  reference  image  W  of  size  KxL  is  chosen  from  the  center  of 

the  field  of  view  of  the  preprocessed  high  resolution  image.  The 
correlation  between  row  i  and  row  i+1  is  given  by 


l  W(i  ,j)  W(i+1  ,j) 


rR(i )  = 


'L  2  1 
T.  W2(i,j)! 

1/2 

f  2 

t.  ir(i+l,j) 

LM  J 

LH  J 

1/2 


(?-12) 


for  i  =  1 ,  2,  3,  . . . ,  K-l  . 
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Table  5. 


Scene 


Target  #1 
Target  H 2 

Position  A 

Target  #3 
Target  #4 


Target  it  1 
Target  #2 

Position  B 

Target  #3 
Target  #4 


Target  #1 
Target  #2 

Position  C 

Target  #3 
Target  #4 


Target  #1 
Target  it 2 

Position  D 

Target  #3 
Target  #4 


Image  registration  using  correlation 
of  adjacent  pixels. 


Correlation  of  adjacent 
pixels  in  each  row  and 
.each  column 

Correlation  of  adjacent 
rows  and  columns 

32x32  Ref. 

32x64  Ref. 

32x32  Ref. 

32x64  Ref. 

X 

- 

3 

- 

X 

X 

X 

1 

X 

X 

X 

2 

X 

X 

1 

1 

X 

X 

4 

X 

X 

X 

X 

X 

X 

X 

X 

1 

i 

c 

X 

2 

4 

X 

1 

1 

X 

4 

X 

X 

2 

1 

1 

1 

1 

1 

1 

2 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

4 

X 

X 

X 

X 
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Step  3: 


Step  4: 


Step  5: 

A  total  or  sixteen  pairs  of  scenes  are  simulated  using  the  above  method 
for  reference  images  of  32  x  32  with  results  presented  in  Table  5.  The 
registration  point  appeared  as  the  first  maximum  for  four  scenes,  and  was 
within  the  first  four  maxima  for  two  additional  scenes.  Registration 
point  was  not  within  first  four  maxima  for  ten  remaining  scenes.  When 
the  reference  image  size  was  increased  to  32  x  64,  the  registration  point 
appeared  as  the  first  maximum  in  the  correlation  surface  for  five  scenes, 
was  within  the  first  four  maxima  for  four  additional  scenes,  and  was  not 
within  the  first  four  maxima  for  the  remaining  six  scenes.  From  Table  5, 
it  is  clear  that  the  performance  of  the  method  based  on  the  correlation 
of  adjacent  rows  and  columns  is  better  than  the  method  based  on  correla¬ 
tion  of  adjacent  pixels  in  each  row  and  each  column.  The  above  method  is 


Similarly,  let  CC(j)  be  the  normalized  correlation  between  col¬ 
umn  j  and  column  j+1  as  given  by 


CC(j) 


K 

l  W(i,j)  W(i , j+1 ) 
i=l 


1/2 

L  2 

7  W  ( 1 , j+1  ) 

1/2 

b=i  J 

J=1  J 

=  12  3 

1  J  4.J  •  •  •  > 

L-l  . 

(2-13) 


The  (K+L-2)-dimensional  feature  vector,  V^,  of  the  digital  image 
W  is  given  by 


\l'u  =  [CR ( 1  )  CR(2)  ....  CR(K-l)  CC(1)  CC(2)  ...  CC(L-l)]  (2-14) 


Feature  vectors  are  computed  for  all  subirnages  of  the  search 
area  and  matched  with  those  of  the  reference. 


27 


expected  to  do  better  if  the  scenes  consist  of  prominent  horizontal  and 
vertical  edges.  A  summary  of  the  results  in  Tables  4  and  5  is  presented 
in  Table  6. 

E.  Image  Registration  Using  Moments  Method 
A  method  of  accomplishing  digital  image  registration  based  on  moments 
of  the  reference  W  and  subimages  of  the  search  area  S  is  given  in  Refer¬ 
ence  1.  For  multiple  image  registration,  it  has  been  shown  that  the 
moments  method  requires  fewer  arithmetic  operations  than  the  standard 
correlation  algorithm.  However,  not  much  is  known  about  the  accuracy  and 
reliability  of  the  method  which  may  be  very  sensitive  to  noise.  The  re¬ 
liability  of  the  method  is  investigated  through  simulation  using  the  same 
16  typical  military  type  scenes  used  in  earlier  simulations.  An  attempt 
is  made  to  accomplish  digital  image  registration  without  intensity  level 
normalization  of  the  window  and  the  search  area  as  described  in  the  steps 
below. 

Step  1:  A  reference  image  of  size  KxL  (32  x  32)  is  taken  from  the  reduced 
high  resolution  FLIR  image  to  include  the  prominent  features  of 
the  image. 

Step  2:  All  moments  of  W  of  order  three  and  less  are  computed  using  the 
relation, 

w  K  L 

m”  =  Z  Y,  rmq  W(£,m)  ,  (2-15) 

pq  £=1  m=l 

for  p,  q  =  0,  1,  2  and  3  such  that  p  +  q  3. 


Table  6.  Summary  of  results  in  Table  4  and  Table 
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Step  3:  All  moments  of  order  3  and  less  are  computed  for  each  subimage 
of  the  search  area.  Let  {tn^}  be  the  moment  sequence  of  the 
subimage  S.  .. 

Step  4:  The  moments  of  the  subimages  of  S  are  matched  with  those  of  the 
reference  image.  The  degree  of  dissimilarity  between  W  and  S.  . 

1  jJ 

is  computed  as 


E(i  )  =  £  (m  w  - 
K  '  '  pq  pq ' 

p,q 


(2-16) 


for  1  <  i  <  M-K+l  ,  1  <  j  <  N-L+l  , 


where  M  and  N  are  the  number  of  rows  and  columns  in  the  search 
area.  Minimum  error  is  expected  at  the  registration  point. 

A  total  of  sixteen  pairs  of  scenes  were  simulated  and  the  method 
failed  for  each  pair.  Even  though  the  reference  images  were  chosen  to 
include  the  prominent  features  of  the  reduced  high  resolution  intayes,  trie 
true  registration  points  did  not  appear  within  the  first  four  minima  in 
the  error  surface.  The  simulation  was  repeated  using  intensity  level 
normalized  windows  and  search  areas  without  any  success.  Therefore,  it 
is  concluded  that  moments  method  is  not  suitable  for  image  registration 
for  the  kind  of  images  considered  in  this  report,  at  least  not  without 
additional  preprocessing  such  as  smoothing  which  was  not  tried. 


F.  Image  Registration  Using  Intraset  and 
Interset  Distances 

A  method  of  accomplishing  multiple  image  registration  using  charac¬ 
teristic  features,  namely  the  intraset  and  interset  distances  is  discussed 
in  detail  in  [!].  Since  features  are  first  extracted  and  matching  is 
performed  on  the  feature  set  only,  this  method  is  computationally  more 
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efficient  than  the  correlation  methods  for  multiple  image  registration. 
However,  the  accuracy  and  reliability  of  feature  matching  algorithms 
depend  very  much  on  the  efficient  selection  of  features.  The  registra¬ 
tion  accuracy  and  reliability  of  this  method  for  bilevel  and  trilevel 
images  are  investigated  through  simulation  using  typical  military  scenes 
described  in  a  previous  section.  The  registration  technique  used  and 
the  simulation  results  are  presented  below. 

Step  1:  A  reference  image  W  of  size  KxL  is  selected  from  the  reduced 
high  resolution  FLIR  image. 

Step  2:  The  mean  n  and  the  standard  deviation  ow  of  the  pixel  values  of 
W  are  computed.  The  window  W  is  segmented  into  three  sets  Sg, 

S^ ,  a  nd  . 

Sg  ,  if  W(a,m)  <  y^|  -  Soy 

W(£,m)  e  Sj  ,  if  W(t5m)  >  +  Sow  (2-17) 

,  otherwise 

for  and  1  £  m  <_  L  . 

Where  S  is  a  scale  factor  greater  than  zero.  The  reference  image 
W  can  be  segmented  into  only  two  sets  SQ  and  S^  by  setting  S 
equal  to  zero. 

Step  3:  The  intraset  distances  of  sets  SQ,  S^ ,  and  S,,  and  the  interset 
distances  between  Sg  and  S^ ,  Sg  and  S,,,  and  S,  and  S£  are  com¬ 
puted  [I].  Let  D  (S.)  and  D  (S . ) ,  i  =  0,  1 ,  2  be  the  x  and  y 
x  i  y  i 

components  of  the  intraset  distances,  and  D  (S.,  S.)  and 

xi  j 

S. ) ,  i  f  j  for  i , .i  =  0,  1,  2  be  x  and  y  components  of 
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the  interset  distances.  The  12  components  of  intraset  and  inter¬ 
set  distances  form  a  12-dimensional  feature  vector  V,,  of  the 
reference  image  W. 

Step  4:  A  simil.  r  feature  vector  V.  -  for  each  subimage  S.  .  of  size  KxL 
of  the  reduced  low  resolution  FL1R  image  (search  area)  is  computed 
according  to  steps  2  and  3. 

Step  5:  The  error  surface  E  is  generated  by  computing  the  Euclidean  dis¬ 
tance  between  V,,  and  V.  .  for  all  subimages  of  the  search  area 
W  i  ,j 

S. 

Step  6:  The  registration  point  is  expected  to  have  minimum  error. 

Reduced  low  resolution  FLIR  images  of  size  406  x  313  and  reduced  high 
resolution  FLIR  images  of  size  98  x  78  described  earlier  in  the  report 
are  used  in  the  : imulation.  For  each  pair  of  scenes  (search  area  and  win¬ 
dow)  the  follow  -g  three  cases  are  considered. 

Case  1:  Refererc  arrays  of  size  16  x  16  are  chosen  to  include  the  prom¬ 
inent  1:  tures  in  the  reduced  high  resolution  FLIR  images.  The 
windows  nd  subimages  of  the  search  area  are  segmented  into  two 
levels  t'  setting  S  equal  to  zero  in  Equation  (2-17).  The  fea¬ 
ture  vec  or,  therefore,  consists  of  one  interset  and  two  intraset 
di stance. . 

Case  2:  Reference  arrays  of  size  24  x  24  are  chosen  to  include  the  p  •">- 
irient  features  in  the  reduced  high  resolution  FLIR  images.  The 
segmentation  of  images  is  the  same  as  Case  1. 

Case  3:  Reference  arrays  of  size  24  x  24  are  chosen  to  include  the  prom¬ 
inent  features  in  the  reduced  high  resolution  FLIR  images  and 
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are  segmented  into  three  levels  Sg,  ,  and  by  setting  S 
equal  to  0. 1  . 

Case  4:  24  x  24  reference  arrays  and  subimages  of  the  search  area  are 

segmented  into  three  levels  Sg,  ,  and  by  setting  S  equal  to 

0.1.  However,  the  set  S£  is  not  considered  in  the  registration 

process.  In  other  words,  the  feature  vectors  consist  only  of 

the  two  intraset  distances  of  sets  Sg  and  ,  and  an  interset 

distance  between  S,.,  and  S, . 

u  i 

The  simulation  results  for  each  of  the  above  cases  is  given  in  detail 
in  the  monthly  reports  of  December  1980  and  January  1981  and  are  summa¬ 
rized  in  Tables  7  and  8.  In  Table  7,  the  appearance  of  the  registration 
point  as  1st,  2nd,  3rd  and  4th  minima  in  the  error  surface  is  indicated 
by  1 ,  2,  3,  and  4,  respectively.  "X"  indicates  that  the  registration 
point  did  not  appear  within  first  four  minima.  A  few  simulations  wt-re 
also  run  for  reference  arrays  of  size  32  x  32  and  32  x  64  with  extremely 
poor  results.  From  this  result  and  the  results  summarized  in  Tables  7 
and  8  it  is  concluded  that  this  method  is  not  accurate  for  the  type  of 
imagery  considered  and  will  not  be  discussed  further. 

G.  A  New  Feature  Matching  Image  Registration  Technique 
A  new  method  of  accomplishing  image  registration  using  features 
based  on  adjacent  pixel  differences  along  the  rows  and  columns  of  a  digi¬ 
tal  image  is  presented.  In  the  feature  matching  technique,  a  digital 
image  is  represented  by  R  features  or  an  R-dimensional  feature  vector, 
and  the  decision-making  process  is  based  on  a  similarity  measure  which, 
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Table  7.  Simulation  results  of  intraset  and 
interset  distance  methods. 


Scene 


Case  1  Case  2  Case  3  Case  4 


Position  A 


Target  #1 
Target  #2 
Target  #3 
Target  #4 


X 

1 

X 

X 


2 

2 

X 

X 


4 

1 

X 

X 


Position  B 


Target  #1 
Target  #2 
Target  43 
Target  44 


X 

4 

X 

X 


X 

4 

3 

X 


X 

X 

1 

X 


X 

2 

4 

X 


Position  C 


Target  41 
Target  42 
Target  #3 
Target  44 


X 

X 

X 

1 


X 

X 

4 

2 


Position  D 


Target  #1 
Target  42 
Target  #3 
Target  44 


1 

X 

X 

X 


X 

2 

X 

X 


4 

3 

2 

X 


X 

3 

2 

X 
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Table  8.  Summary  of  intraset  and  interset 
simulation  results  presented  in  Table  7. 


Registration  Point  is  Registration  Point  is  Registration  Point 
the  first  minimum  2nd,  3rd,  or  4th  not  in  first 

Minimum  four  minima 

Case  16  1  9 

Case  22  3  11 

Case  31  7  8 

Case  43  6  7 
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in  turn,  is  expressed  in  terms  of  a  distance  measure  (Euclidean  distance) 
or  the  correlation  of  feature  vectors  (cosine  of  the  angle  between  the 
two  vectors). 

Let  W  be  a  digital  image  (window)  of  size  KxL  whose  pixels  can  as¬ 
sume  one  of  G  possible  levels  on  the  gray  scale 

0  <  W(i,j)  <  G-l 

for  1  _<  i  £  K  and  1  <  j  <  L  . 

The  feature  elements  of  each  row  and  column  of  W  are  computed  using  equa¬ 
tions  (2-18)  and  (2-19),  respectively. 


z  [W(k, j )  -  W(k,j+1 )] 

^  - - - ,  for  k  =  1,  2,  ....  K  (2-18) 

L  W(k,j)2 


t  [W ( i  ,£.)  -  W(i+1  ,st)]2 

C£  =  ~ - K - »  for  1  =  2>  •••>  L  (2'19) 

l  W(i,fc)2 
i=l 

The  feature  element  R.  (C  )  is  always  greater  than  or  equal  to  zero,  as- 

K 

suming  zero  value  only  when  all  the  pixels  in  the  kth  row  (nth  column)  of 
W  are  at  the  same  level  on  the  gray  scale. 

A  computationally  more  efficient  way  of  extracting  row  and  column 
features  of  a  digital  image  of  size  KxL  is  to  use  Equations  (2-20)  and 
(2-21)  instead  of  (2-18)  and  (2-19),  respectively. 


Rk  -  1=L 


E  |W(k,j)  -  W(k,j+l)| 


,  for  k  =  1  ,  2,  .  , . ,  .K 


(2-20) 


E  |W(k,j)| 
j=l 


E  |W\i  ,P.)  -  W(i+l,a)| 

c  _ _ 

S  K 


,  for  z  =  1,  2,  ....  L  (2-21) 


r  j W( i ,?.)  | 


The  ( K+L ) -d i mensi ona T  feature  vector,  V,,  of  the  digital  image  W  is  now 
formed  as  given  by 


»I  ’  R2  •••  «KC1  C2  V 


(2-22) 


This  feature  vector  is  correlated  with  the  feature  vectors  of  the  sub¬ 


images  of  size  KxL  of  the  search  area  to  determine  the  registration  point. 


1.  Derivation  of  the  Computational  requirements 


for  Multiple  image  registration 


The  number  of  additions  and  multiplications  required  to  implement 


the  new  feature  matching  algorithm  for  a  search  area  of  size  MxN  and  n 


reference  images  (windows)  of  size  KxL  are  derived  in  the  following  steps: 


1,  Computation  of  R.  as  given  by  Equation  (2-20)  requires  (3L-4)  additions 

K 

(subtractions)  and  one  multiplication  (division).  Similarly,  (3K-4) 


additions  and  one  multiplication  are  needed  to  compute  Cfc.  Since 


there  are  K  rows  and  L  columns,  the  total  number  of  additions  and 


n 


IH 
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multiplications  required  to  compute  the  feature  vector  of  a  digi¬ 
tal  image  W  of  size  KxL  are  given  by 

Total  number  of  additions  =  K(3L-4)  +  L(3K-4) 

=  6KL  -  4 ( K+L  ) 

Total  number  of  multiplications  =  K+L 

2.  In  order  to  compute  the  feature  vector  V.  .of  the  subimage  S.  .also 

i,j  i,j 

of  size  KxL,  6KL  -  4(K+L)  additions  and  (K+L)  multiplications  are  re¬ 
quired  as  outlined  in  step  1. 

3.  The  feature  vector  V-  ~  of  the  subimage  S.  9  can  be  computed  with  very 

few  arithmetic  operations  once  V.  , ,  the  feature  vector  of  subimage 

i ,  * 

S.  , ,  is  computed. 

1  *  » 

Now, 


L-l 


Rkd,l) 


T  |si  ,(k,j)  -  S  (k,j+l)| 

1=1  '  >  1  I »  I 


ill 


i  |S.  (k,j)| 

j=l  1,1 


Nk(i,l) 

DjTTrrr 


(2-23) 


for  k  =  1 ,  2,  ....  K  . 


where 


1.-1 

Nk(i,l)  =  I 
k  j=l 


|Sifl(k,j)  -  Sisl(k,j+1 ) |  ,  and 


L 

D.(i,l)  =  z  |S  1  (k,  j )  [ 

k  j=1  I,' 
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Then 


Rk(i,2) 


L-l 

z  Is.  -  s.  2(k,j+i)| 

j=1  _ _ 111 _ 

I. 

e  ■  is.  2(kj}| 

.M  11 


Nk(i,l)  -  |Su(k,l)-S.  ^  1  ( k ,  2 )  |  +  |Si#2(k,M)-Si>z(K,L) 
Ok(i,D  -  I  si  -j  (k,  T  )  |  +  [S.  2(k,L)| 

(2-24) 


It  is  seen  from  Equation  (2-24),  that  R ^ ( i , 2 )  can  be  computed  from 
Rk(i,l)  and  in  general  Rk(i,j)  can  be  computed  from  Rk(i,j-i)  by  per¬ 
forming  only  6  equivalent  addition  and  one  multiplication  operations. 
Computation  of  the  K  row-feature  elements  of  5.  0,  therefore,  requires 
6K  additions  and  K  multipMcations.  Also,  since 

Cj, (i «2)  =  C£+1  (i ,1 )  ,  for  £  -  1,  2 . L-l  , 

it  is  only  necessary  to  compute  Ck(i,2)  which  requires  (3K-4)  addi¬ 
tions  and  one  multiplication  as  shown  in  Step  1.  The  computation  of 
the  feature  vector  V.  ?  from  V.  therefore,  requires  a  total  of 
(9K.-4)  additions  and  (K+l)  multiplications. 

4.  There  are  (M-K+l ) (N-L+l )  reference  points  or  subimages  in  a  search 
area  of  size  MxN.  The  computation  of  the  feature  vector  associated 
with  the  reference  point  (1,1)  requires  6KL  -  4(K+L)  additions  and 
K+L  multiplications  (Step  1).  The  computation  of  the  remaining  (N-L) 
reference  points  in  the  first  row  requires  a  total  of  (9K-4)(N-L)  ad¬ 
ditions  and  (K+l) (N-L)  multiplications.  The  computation  of  the  fea¬ 
ture  vectors  of  all  the  (M-K+l ) (N-L+l )  subimages,  therefore,  requires 
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(M-K+l  )[6KL-4(K+L)  +  ( 9K-4 ) (N-L ) ]  additions,  and 

(M-K-fl  )[(K+L)  +  ( K+l ) (N-L ) ]  multiplications  .  (2-25) 

5.  The  correlation  surface  is  generated  by  correlating  the  reference 
feature  vector  with  each  of  the  (M-K+l ) (N-L+l )  feature  vectors  of 
the  subimages  using  the  relation 


Ci » J 


iv 


for  all  i , j 


W1 


|vij 


(2-26) 


Computation  of  C.  .  for  feature  vectors  of  dimension  (K+L)  requires 
i  >  J 

(3K+3L-3)  additions  and  (3K+3L+5)  multiplications.  The  generation  of 
the  correlation  surface,  therefore,  requires  a  total  of 

(M-K+l ) (N-L+l ) (3K+3L-3)  additions  and 
(M-K+l  ) (N-L+l ) ( 3K+3L+5)  mul tipi i cations. 

6.  Assuming  that  there  are  n  windows,  the  multiple  image  registration 
process  requires  Steps  1  and  5  to  be  performed  n  times.  Therefore, 
for  multiple  image  registration. 

Total  number  of  additions  required  = 

-  n[6KL-4(K+L)] 

+  (M-K+l )[6KL-4(K+L)+(9K-4) (N-L)] 

+  n (M-K+l ) (N-L+l ) (3K+5L-3) 

Total  number  of  multiplications  required  = 

=  n ( K+L ) 

+  (M-K+l  ) [(K+L)+(K+1 ) (N-L) ] 

+  n (M-K+l ) (N-L+l ) (3K+3L+5)  (2-27) 
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for  M  -  240,  N  =  256,  K  =  32,  arid  L  =  32 
Total  number  of  additions  =  8893613  n  +  14526336 
Total  number  of  multiplications  =  9263989  n  +  1558304 
The  number  of  additions  and  multiplications  which  are  required  to  im¬ 
plement  correlation,  SSDA,  moments  method  and  the  new  feature  match¬ 
ing  method  for  various  values  of  n,  are  shown  in  Figures  (2-5)  and 
(2-6),  respectively.  Assuming  that  multipl ication,  division,  and 
squaring  operations  are  equivalent  and  that  the  time  required  to  per¬ 
form  one  real  multiplication  is  three  times  that  required  to  perform 
one  real  addition,  the  total  number  of  equivalent  real  additions  re¬ 
quired  to  implement  the  various  methods  are  shown  in  Figure  [1-1). 

It  can  be  concluded  from  Figures  <2-5)  through  (2-7)  that,  the  new 
feature  matching  technique  is  computationally  more  efficient  than 
any  of  the  template  matching  or  feature  matching  methods  considered 
so  far  for  multiple  image  registration. 

2.  Simulation  Results 

An  effort  is  made  to  study  the  performance  and  registration  accuracy 
of  the  new  feature  matching  method  through  simulation.  Reduced  low  reso¬ 
lution  FUR  images  (search  area)  of  size  406  x  313  and  reduced  high  reso¬ 
lution  FUR  images  of  size  98  x  78  are  used  for  simulation.  The  registra¬ 
tion  technique  and  the  simulation  results  are  presented  below  for  all  the 
cases  considered. 


16.8x70 


Figure  2-5.  Number  of  additions  to  register  n  windows. 


Moments 


Number  of  equivalent  additions  to  register  n  windows 
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Case  1:  Reduced  low  and  high  resolution  FLIR  images  a^e  used  for  the 

simulation  purpose  without  intensity  level  normalization  or  any 
other  preprocessing.  The  registration  technique  is  described 
below. 


.Step  1:  A  reference  image  W  of  size  KxL  is  selected  from  the  reduced 
high  resolution  FLIR  image. 

Step  2:  The  feature  elements  of  each  row  and  column  of  W  are  computed 
using  Equations  (2-18)  and  (2-19),  respectively. 

Step  3:  The  (K+L)-dimensional  feature  vector,  Vw>  of  the  reference  image 
W  is  now  formed  as  given  by 


Step  4:  A  subarray  of  size  120  x  240  which  includes  the  target  area  of 

interest  is  extracted  from  the  406  x  313  reduced  low  resolution 

FLIR  image  and  is  used  as  the  search  area. 

Step  5:  The  (K+L)-dimensional  feature  vectors,  V.  for  each  subimage 

i  >  J 

S.  .  of  size  KxL  is  computed  as  in  Steps  2  and  3  and  the  cosine 

1  9  J 

of  the  angle  between  the  two  feature  vectors  V,.  and  V.  .is  com- 

w  i  »J 

puted  using  the  relation 


for  all  i,j 


Step  6:  C-  .  is  expected  to  have  a  maximum  at  the  registration  point. 

1  1  J 

Simulations  were  run  for  a  reference  array  of  size  32  x  32, 
choosing  the  reference  from  the  center  of  the  field  of  view  of 
the  reduced  high  resolution  FLIR  image. 
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Case  2:  As  in  Case  1  with  32  x  64  reference  images. 

Case  3:  Case  3  is  similar  to  Case  1  except  that  the  64-dimensional  fea¬ 
ture  vector,  Vw,  of  the  reference  image  W  of  size  32  x  32  is 
formed  by  computing  the  row  and  column  feature  elements  using 
Equations  (2-20)  and  (2- 21)  instead  of  Equations  (2-18)  and 
(2-19),  respectively.  The  feature  vectors  V.  .  for  each  sub- 

image  S.  .  of  size  32  x  32  are  also  computed  using  Equations 
■  >J 

(2-20)  through  (2-22). 

Case  4:  As  in  Case  3  with  32  x  64  reference  images. 

Case  5:  Since  the  two  images,  the  window  W  and  its  matching  subimage 

S.  .of  the  search  area  are  obtained  from  different  sensors,  the 

1  5J 

corresponding  pixels  in  the  two  images  may  not  have  the  same 

pixel  value  due  to  the  difference  in  gain  and  dc  bias  of  the 

two  imaging  systems.  In  such  a  case  it  is  reasonable  to  assume 

that  S.  .  and  W  are  related  by  the  following  equation 
I  *J 

W(s,,m)  =  a  5.  . (*.,m)  +  b  ,  1  £  1  <  K  and  1  <  m  <  L 

where  a  (gain)  and  b  (dc  bias)  are  constants.  Since  the  numera¬ 
tors  of  Equations  (2-20)  and  (2-21)  are  computed  by  taking  the 
difference  between  adjacent  pixels,  they  are  independent  of  the 
constant  b.  Since  the  denominators  are  functions  of  the  constant 

b,  the  matching  of  feature  vectors  of  W  and  S.  .  may  lead  to 

9  J 

false  r  nstration.  This  false  registration  can  be  avoided  by 
preprocessing  the  two  images  so  as  to  have  zero  means  or  equal 
means.  The  constant  multiplication  factor  "a"  does  riot  result 
in  false  registration  because  the  decision  making  process  is 
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based  on  a  similarity  measure  which  is  expressed  as  the  correla¬ 
tion  of  feature  vectors  (see  the  Equation  in  Step  5). 

The  registration  process  is  identical  to  that  of  Case  3,  ex¬ 
cept  that  both  the  high  resolution  and  low  resolution  FLIR  images 
are  preprocessed  using  a  moving  area-average  filter  before  the 
extraction  of  the  64-dimensional  feature  vector.  Preprocessing 
of  the  two  images  is  done  using  an  algorithm  which  replaces  the 
value  of  each  pixel  by  its  value  minus  the  average  pixel  value 
of  the  11  x  11  array  centered  about  the  pixel.  The  two  prepro¬ 
cessed  images  will  then  have  zero  or  nearly  equal  means.  Simula¬ 
tion  results  indicated  a  significant  improvement  in  the  perfor¬ 
mance  when  preprocessed  images  are  used. 

Case  6:  Reference  image  of  size  32  x  64  chosen  from  the  center  of  the 
field  of  view  of  the  reduced  high  resolution  FLIR  image  is  used 
for  the  simulation  purpose.  The  registration  technique  is 
similar  to  Case  3.  Both  the  high  resolution  and  low  resolution 
FLIR  images  are  preprocessed  using  a  moving  area-average  filter 
described  in  Case  5  before  the  extraction  of  the  9G-dimensicnal 
feature  vector. 

Case  7:  The  registration  technique  is  similar  to  Case  4  except  that  the 
feature  vector  is  computed  in  a  slightly  different  way.  The 
32  x  64  array  is  processed  as  two  separate  blocks  of  size  32  x 
32  and  the  registration  process  is  describee  below. 

Step  1:  The  reference  image  W  of  size  32  x  64  is  divided  into  two  sub- 
images  and  of  size  32  x  32  each. 
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Step  2:  The  feature  elements  of  each  row  and  column  of  and  are 
computed  using  Equations  (2-28)  and  (2-29),  respectively. 


32-1 

T  i [W, (k, j )  -  W  (k, j+1 )] | 

R  _  i=J _ _ 

32  (2-28) 

2  I W. (k,  J ) I 
j=l  £ 


32-1 

l  1 [W  (i ,k)  -  W  (i+1 ,k)][ 

r  =  i=1  _ _ _ 

Lka  32 

£  |W  (i,k)| 
i=l 

for  k  =  1,  2,  3,  32,  and  £=1,2. 


(2-29) 


Step  3:  The  128-dimensional  feature  vector,  V^,  of  the  reference  image 
W  is  now  formed  as  given  by 


VW  =  ER1,1  R2,l  R32,l  R1 ,2  R2,2  ”•  R32,2  C1 ,1  C2,l  •"  C32,l 
C1 ,2  C2,2  •••  C32,2^ 


Step  4:  A  128-dimensional  feature  vector,  V.  .,  for  each  subimage  of 

1  5  J 

size  32  x  64  from  the  search  area  is  computed  according  to  Steps 
2  and  3. 

Step  5:  The  two  feature  vectors  V..  and  V.  .  are  correlated  to  find  the 

w  1  ,j 

registration  point. 

The  simulation  results  for  each  of  the  above  cases  is  given  in  de¬ 
tail  in  the  monthly  reports  of  May  through  July  1981  and  are  summarized 
in  Tables  9  and  10.  The  simulation  results  indicate  that  the  new  feature 
matching  technique  performs  better  than  the  binary  correlation  method 
and  at  the  same  time  requires  less  computation  since  the  feature  vectors 
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Table  9.  Simulation  results  of  new  feature  matching  image 


registration 

technique. 
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Table  10.  Summary  of  results  in  Table 
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are  generated  recursively.  The  performance  is  still  improved  with  the 
modified  features  where  the  absolute  value  of  the  difference  between  ad¬ 
jacent  pixels  is  taken  instead  of  the  squared  value  resulting  in  further 
reduction  in  computation.  Preprocessing  of  the  images  using  a  moving 
area-average  filter  before  the  extraction  of  the  features  resulted  in  a 
significant  improvement  in  the  registration  accuracy  of  the  modified  fea¬ 
ture  matching  technique.  Further  improvement  in  the  registration  accuracy 
when  a  32  x  64  image  array  is  processed  as  two  separate  blocks  of  size 
32  x  32  is  evident  from  the  results  summarized  in  Table  9.  Based  on  the 
above  facts  it  can  be  concluded  that  for  the  type  of  military  scenes  con¬ 
sidered  the  new  modified-feature  matching  technique  with  preprocessed 
images  performs  better  than  all  the  methods  considered  so  far,  where  the 
reference  is  taken  from  the  center  of  the  FOV,  and  at  the  same  time  re¬ 
quires  less  computation  time,  thus  making  the  method  more  attractive  for 
real-time  implementation. 

H.  Image  Registration  Using  Haar  Coefficients 
Orthogonal  transforms,  such  as  the  Haar  transform  (HT),  convert  a 

i 

digital  image  in  the  "spatial -domain"  into  a  "transform-domain"  image 
that  contains  the  same  information  in  a  nonspatial  representation.  The 
Haar  transform  is  derived  from  the  Haar  matrix  [6]  whose  elements  are 
either  1,  -1,  or  0  multiplied  by  powers  of  J2.  An  8  x  8  Haar  matrix  is 
shown  in  Figure  2-8. 


\ 

\ 

\ 
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r  iJT, 
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Figure  2-8.  An  8  x  8  Haar  matrix. 


The  two  dimensional  Haar  transform  [F]  of  a  digital  image  [f]  is  defined 
by 

[F]  =  H[f]  HT  (2-30) 

where  [F],  [f],  and  H  are  all  of  the  same  order. 

The  Haar  transform  can  be  visualized  as  a  course-to-f ine  sampling 
process  of  the  input  data  sequence,  where  the  first  transform-coefficient 
represents  the  mean  value  of  the  components  of  the  sequence,  the  second 
represents  the  average  difference  of  the  first  N/2  components  and  the 
second  N/2  components,  and  so  on,  with  the  last  N/2  transform-coefficients 
measuring  the  adjacent  differences  of  the  data  elements  taken  two  at  a 
time.  The  Haar  transform  coefficients  may  be  called  the  sequency  (spatial 
frequency)  components  and  are  analogous  to  the  Fourier  transform  coeffi¬ 
cients  which  are  called  the  frequency  components.  The  sequency  of  any 
row  (or  column)  of  a  Haar  matrix  is  the  number  of  sign  changes  along  that 
row  (or  column)  [10].  One  of  the  several  properties  of  Haar  transform  is 
its  ability  to  store  a  large  percentage  of  the  image  information  in  a  few 
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coefficients,  so  that  a  feature  vector  with  only  a  small  number  of  ele¬ 
ments  is  required  for  retention  of  significant  image  information.  The 
criterion  for  the  selection  of  features  for  matching  purposes  may  be  one 
of  the  following: 

1.  Select  all  transform  coefficients  whose  absolute  values  are 
greater  than  a  predetermined  threshold. 

2.  Select  those  transform  coefficients  which  have  the  largest 
variation  or  variance  across  image  classes. 

3.  Select  a  predetermined  number  of  highest  transform  coefficients. 

4.  Choose  the  low  sequency  components  assuming  that  a  major  por¬ 
tion  of  the  high  sequency  components  represent  noise. 

Computatio,!  of  the  two  dimensional  Haar  transform  of  a  digital  image 
can  be  best  illustrated  by  using  the  two  dimensional  basis  plane  as  shown 
in  Figure  2-9  for  an  8  x  8  transform.  Each  of  the  transform  coefficients 
may  then  be  considered  proportional  to  the  correlation  of  the  image  and 
the  corresponding  basis  plane.  For  example,  the  (3,  4)  coefficient  in 
the  transform-domain  image  is  computed  as  follows.  Every  element  of  the 
8x8  spatial -domain  image  of  Figure  2-10  is  multiplied  by  the  correspond¬ 
ing  element  in  the  (3,  4)  basis  plane  shown  in  Figure  2-11  and  the  pro¬ 
ducts  are  added  together.  The  resulting  product  is  then  multiplied  by 
the  proper  weighting  factor  which  is  2  for  the  (3,  4)  coefficient  as  in¬ 
dicated  in  Figure  2-12.  The  (3,  4)  coefficient  of  the  transform  domain 
image  is  therefore  given  by 

F(3,  4}  =  2(f 1 5  +  f16  “  f]7  -  f18  +  1 25  +  f26  "  f27  ~  f28 
'  f35  '  f  36  +  f37  +  f38  -  f45  “  f46  +  f47  +  W  ’ 


Figure  2-9.  Basis  Plane  for  an  8x8  Haar  transform. 
B'lack=+1;  crosshatched=-l ;  white=0. 


56 


1 

1 

4i 

ft 

2 

2 

2 

2 

■ 

■ 

2 

2 

2 

2 

- 

K 

D 

2 

2 

2/2 

2/2 

2/2 

2/2 

! 

72 

J2 

2 

2 

2/2 

zft 

2/2 

_ 

2/2 

2 

2 

2/2 

2/2 

■ 

B 

B 

B 

2 

2 

2/2 

2/2 

■ 

B 

B 

B 

2 

2 

2/2 

zjz 

B 

B 

fl 

1 

2 

2 

2/2 

zft 

B 

B 

B 

B 

I 

x 


Figure  2-12.  The  weighting  factors  for  the  8x8 
transform-domain  image. 
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Other  coefficients  are  computed  in  a  similar  way  by  repeating  the  process 
of  correlation  with  the  corresponding  basis  plane. 

The  Haar  transform  can  he  computed  efficiently  by  using  a  fast  Haar 
transform  (PUT).  A  flow  diagram  for  computing  an  8-point  fast  Haar  trans¬ 
form  is  shown  in  Figure  2-13.  In  the  diagram  the  input  sequence  numbers 
are  represented  by  t..  and  the  transformed  coefficients  by  F...  Each  node 
represents  a  summing  junction  and  multiplication  by  constants  at  various 
stages  of  the  computation  is  indicated  by  -1  or  /2  or  2  along  the  branches. 
A  two  dimensional  Haar  transform  is  then  computed  by  repeated  application 
of  the  one  dimensional  FHT. 

1 .  Modified  and  Rationalized  Haar  Transforms 

Since  some  of  the  elements  in  the  Haar  matrix  are  irrational  numbers 
(powers  of  /2),  a  variation  of  the  transform  called  the  modified  Haar 
transform  (MHT),  which  does  not  require  a  multiplier,  is  more  attractive 
for  digital  implementation.  A  modified  Haar  matrix  whose  elements  are 
either  1,  -1,  or  0  is  shown  in  Figure  2-14. 

Another  variation  of  the  Haar  transform  that  does  not  involve  num¬ 
bers  which  are  powers  of  «/2  is  the  rationalized  Haar  transform  (RHT) .  A 
rationalized  Haar  matrix  has  elements  which  are  either  1,  -1,  0  or  inte¬ 
ger  powers  of  2  as  shown  in  Figure  2-15.  The  multiplication  by  pow^rs- 
of-two  is  easily  implemented  as  binary  shifts. 

All  the  transform  techniques  described  possess  a  fast  computational 
algorithm  when  the  order  of  the  transform  being  computed  is  an  integer 
power-of-two.  The  computational  requirements  of  the  FHT,  MHT,  and  RHT 
are  compared  with  several  other  well  known  transforms  in  Table  11. 
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Figure  2-15.  An  8  x  8  rational i zed-Haar  matrix. 


Both  modified  Haar  and  rationalized  Haar  transforms  possess  the  fol¬ 
lowing  important  properties  of  the  Haer  transform: 

•Images  can  be  reconstructed  by  taking  the  inverse  Haar  transform 

•They  possess  a  fast  computational  algorithm 

•Transform  coefficients  are  locally  as  well  as  globally  sensitive 


utaticnal  requirements  of  various  two-dimensional  transforms. 
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•Transforms  are  capable  of  representing  a  digital  image  with  a  few 
constituent  elements  in  the  transform  domain  to  a  high  degree  of 
accuracy 

In  a  transform  domain  image,  the  low  sequency  large  amplitude  components 
represent  the  image  background  and  the  high  sequency  large  amplitude  com¬ 
ponents  represent  the  image  edges.  Since  the  high  sequency  components  cf 
the  modified  Haar  transform  are  weighted  less  than  the  corresponding  com¬ 
ponents  of  the  Haar  transform,  the  MHT  can  be  considered  as  a  kind  of 
low-pass  filter  with  a  smoothing  effect.  The  RHT,  because  the  high  se¬ 
quency  components  are  weighted  more,  tends  to  represent  well  defined 
edges  in  the  image  domain  as  edges  in  the  transform  domain  with  large 
component  values.  Therefore,  RHT  produces  an  enhancing  effect  on  the 
image  edges. 

2 .  Selection  of  Haar  Coefficients  as  Features 

One  of  the  several  criteria  for  selection  of  features  for  matching 
purposes  is  to  select  those  transform  coefficients  which  have  the  largest 
variation  or  variance  across  image  classes.  The  computation  of  the  vari¬ 
ance  of  the  coefficients  of  the  transformed  image  of  size  KxK  using  a 
FLIR  image  of  size  MxN  requires  (M-K+l )(N-K+1 )  Haar  transforms  to  be  com¬ 
puted  to  include  all  possible  image  classes.  Since  it  is  extremely  time 
consuming  to  compute  (M-K+l ) (N-K+l )  Haar  transforms  for  1  arge  M  and  N, 
only  a  small  number  of  equally  spaced  subimages  of  size  KxK  are  selected 
from  the  FLIR  image  and  used  for  the  computation  of  the  coefficient  vari¬ 
ances.  The  computation  of  the  transform  coefficient  variances  for  reduced 
lew  resolution  FLIR  images  is  explained  in  the  following  steps: 


1.  Four  hundred  equally  spaced  subimages  of  size  32  x  32  are  selected 
from  the  reduced  low  resolution  FUR  image  of  size  406  x  313. 

2.  The  Haar  transform  of  each  of  the  subimages  selected  in  Step  1  is 
computed. 

3.  Steps  1  and  2  are  repeated  for  four  different  low  resolution  FUR 
images,  thus  yielding  a  total  of  1600  transformed  images. 

4.  Variance  of  each  transform  coefficient  is  then  computed  over  the 
1600  transformed  images  of  Step  3. 

5.  Steps  1  through  4  are  repeated  using  modified  Haar  transform  and 
rationalized  Haar  transform. 

The  variance  distribution  of  the  50  transform  coefficients  having 
the  largest  variation  across  image  classes  is  shown  in  Figure  2-16  for 
the  Haar,  modified  Haar  and  rationalized  Haar  transformed  images.  The 
locations  of  the  64  largest  variance  transform  coefficients  for  the  Haar, 
modified  Haar  and  rationalized  Haar  transformed  images  are  shown  in  Fig¬ 
ures  2-17,  2-18  and  2-19,  respectively.  In  these  fiyures,  "1"  indicates 
the  location  of  the  coefficient  with  the  largest  variance,  "2"  indicates 
the  location  of  the  coefficient  with  the  second  largest  variance  and  so 
on.  The  transform  coefficient  variances  for  reduced  high  resolution  FLIR 
images  are  computed  in  a  similar  way.  The  1600  subimages  of  size  32  x  32 
are  chosen  by  selecting  100  equally  spaced  subimages  from  each  of  the  16 
different  reduced  high  resolution  FUR  images  of  size  98  x  78.  The  dis¬ 
tribution  of  the  transform  coefficients  with  the  50  largest,  variances  is 
shown  ir,  Figure  2-20  and  the  locations  of  the  64  largest  variance  trans¬ 
form  coefficients  for  the  Haar,  modified  Haar  and  rationalized  Haar 


Rationalized  rlaar  transform 
Haar  transform 


igure  2-16.  Transform  coefficient  variance  distribution  for  low  resolution  FUR  images. 


Figure  2-1/.  Location  of  the  64  Haar  transform  coefficients  having 
largest  variances  (Low  resolution  FLIR  images). 


Location  of  the  64  rationalized  Haar  transform  coefficients 
haviny  largest  variances  (Low  resolution  FL1R  images). 
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transformed  images  are  shown  in  Figures  2-21,  2-22,  and  2-23,  respectively. 
From  Figures  2-16  through  2-23  the  following  conclusions  are  made: 

].  a  few  transform  coefficients  (about  30)  have  relatively  large 

variances. 

2.  The  low  sequency  components  of  the  transformed  image  have  higher 
variances  compared  to  those  of  the  high  sequency  components  for 
the  scenes  used. 

The  above  two  facts  can  be  used  to  reduce  the  dimensionality  of  the  fea¬ 
ture  vector  while  retaining  the  prominent  features  of  the  image  for 
matching  purposes. 


I 

%. 

'rt- 


2.  To  compute  the  Haar  transform  of  a  window  of  size  KxK,  an  additional 
4K(K-1)  additions  are  needed. 


Figure  2-21.  Location  of  the  64  Haar  transform  coefficients  having 
largest  variances  (High  resolution  FLIR  images). 


Figure  2-22.  Location  of  the  64  modified  Haar  transform  coefficients 
having  largest  variances  (High  resolution  FLIR  images). 


Figure  2-23.  Location  of  the  64  rationalized  Haar  transform  coefficients 
having  largest  variances  (High  resolution  FLIR  images). 
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3.  The  correlation  surface  is  generated  by  correlating  the  reference 
feature  vector  with  each  of  the  (M-K+l ) (N-K+l )  feature  vectors  of 
the  subimage'  using  the  relation 


lv»  VP 


for  all  i ,j 


Vi,j 


Computation  of  C.  for  feature  vectors  of  dimension  P  requires 
(3P-3)  additions  and  (3P+5)  multiplications.  The  generation  of  the 
correlation  surface,  therefore,  requires  a  total  of  (M-K+l ) (N-K+l ) 
(3P-3)  additions  and  (M-K+l ) (N-K+l ) (3P+5)  multiplications. 

4.  Assuming  there  are  n  windows,  the  multiple  ink...  registration  process 
requires  Steps  2  and  3  to  be  performed  n  times.  Therefore, 

Total  number  of  additions  =  4K(K-1 ) (M-K+l ) (N-K+l ) 

+  n[4K(K-l)  +  (3P-3) (M-K+l) (N-K+l)] 

Total  number  of  multiplications  =  n(M-K+l ) (N-K+l ) (3P+5) 

For  M  =  240.  N  =  256,  K  =  32,  and  P  =  64 

TOv.al  number  of  additions  =  186595200  +  8891693  n 

Total  number  of  multiplications  -  9263925  n 

The  number  of  multiplications  required  is  negligible  when  compared 
co  the  number  of  additions.  The  numbe1"  of  additions  and  multiplications 
which  are  required  to  implement  correlation,  SSDA,  moments  method  a nd  Haar 
transform  technique  for  various  values  of  n,  are  shown  in  Figures  2-24  and 
2-25  respectively.  Assuming  that  multiplication,  division  and  squaring 
operations  are  equivalent  and  that  the  time  required  to  perform  one  real 
multiplication  is  three  times  that  required  to  perform  one  real  addition 


14.4x10 


Figure  2-24.  Number  of  additions  to  register  n  windows 


(IBM  370),  the  total  number  of  equivalent  real  additions  required  to  im¬ 
plement  the  various  methods  is  shown  in  Figure  2-26. 

4 .  Simulation  of  Haar  Transforms 

Haar  transform  and  the  variations  of  the  Haar  transform  tend  to  re¬ 
sult  in  uncorrelated  coefficients  in  the  transformed  image  (to  decompose 
the  image  into  a  relatively  small  number  of  large  amplitude  coefficients 
representing  most  of  the  image  information)  of  an  originally  highly  cor¬ 
related  image.  These  uncorrelated  coefficients  are  expected  to  provide 
a  set  of  features  which  are  more  efficient  for  matching  purposes  than  a 
similar  number  of  features  in  the  originally  highly  correlated  digital 
image.  However,  the  performance  and  accuracy  of  the  image  registration 
algorithm  depends  largely  on  the  selection  of  the  transform-coefficients 
for  matching  purposes.  Several  simulations  were  run  to  study  the  effi¬ 
ciency  and  performance  of  the  image  matching  technique  using  Haar  trans¬ 
form-coefficient  feature  vectors.  Reduced  low  and  high  resolution  FLIR 
images  described  earlier  are  used  in  the  simulation.  The  registration 
technique  and  the  simulation  results  are  presented  below  for  all  the 
cases  considered. 

Case  1:  The  registration  technique  is  described  in  the  following  steps. 
Step  1:  A  reference  image  W  of  size  32  x  32  is  selected  from  the  reduced 
high  resolution  FLIR  image  so  as  to  include  prominent,  fea¬ 
tures  of  the  image. 

Step  2:  The  two  dimensional  Haar  transform  of  W  is  computed  by  using  the 
fast  Haar  transform  (FHT). 

Step  3:  A  64-dirnensional  feature  vector  V'w  consisting  of  the  64  highest 

coefficients  of  the  transformed  image  is  formed  after  setting  the 


s  method  (moments  of  order  <  3) 


equivalent  additions  to  register  n  windows 
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Step  4 


Step  5 


Step  6 


Step  7 
Case  2 


(1,1)  element  of  the  transformed  image  to  zero.  Since  the  (1,1) 
element  represents  the  mean  of  the  pixel  values  of  W,  setting 
this  element  to  zero  simply  ignores  the  image  mean  values  which 
may  be  different.  All  other  coefficients  are  insensitive  to  the 
differing  means  since  they  are  the  sums  and  difference  of  pixel 
values  (i.e.,  subtracting  the  mean  value  before  computing  these 
coefficients  will  not  change  the  coefficient  values). 

In  order  to  reduce  simulation  computation  time,  a  subarray  of  size 
100  x  180  which  includes  the  target  area  of  interest  is  extracted 
from  the  406  x  313  reduced  low  resolution  FLIR  image  (search 
area) . 

The  Haar  transform  of  each  subimage  of  size  32  x  32  from  the 
search  area  is  computed  and  a  64 -dimensional  feature  vector,  V-  • ; 

1  so 

is  formed  for  each  subimage  in  such  a  way  that  the  elements  of 
V.  .  and  V.  have  common  locations  in  the  transformed  images. 

1  aJ  W 

The  cosine  of  the  angle  (phase  correlation)  between  the  two  fea¬ 
ture  vectors  Vw  and  V..  ^  is  computed  using  the  relation 


Ci,j 


VT  v  . 
w  Vl,j 


ngi  iivijii 


,  for  all  i , j 


C.  .is  expected  to  have  a  maximum  at  the  registration  point. 

In  Case  1,  the  reference  images  are  chosen  to  include  promi¬ 
nent  features  of  the  image  and  was  not  necessarily  taken  from  the 
center  of  the  field  of  view  cf  the  reduced  high  resolution  FLIR 
image.  Case  2  is  identical  to  Case  1,  except  that  the  reference 
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images  are  chosen  from  the  center  of  the  field  of  view  of  the 
reduced  high  resolution  FLIR  image. 

Case  3:  In  this  case,  reference  arrays  of  size  32  x  64  are  chosen  from 
the  center  of  the  field  of  view  of  the  reduced  high  resolution 
FLIR  image.  Since  the  computation  of  the  FHT  requires  the  image 
array  to  be  square,  the  32  x  64  reference  array  is  processed  as 
two  separate  blocks  of  size  32  x  32  each.  The  registration  al¬ 
gorithm  is  described  below. 

Step  1:  The  selected  reference  image  of  size  32  x  64  is  divided  into 
two  subimages  and  of  size  32  x  32  each. 

Step  2:  Two  dimensional  Haar  transforms  of  VI ^  and  are  computed  sepa¬ 
rately  using  the  FHT. 

Step  3:  A  128-dimensional  feature  vector  Vw  consisting  of  the  64  highest 
coefficients  of  the  transformed  image  of  and  64  highest  coef¬ 
ficients  of  the  transformed  image  of  W2  is  formed.  The  (1,1) 
elements  of  the  transformed  images,  which  represent  the  means 
of  the  pixel  values  of  W1  and  V^,  are  not  considered  in  the 
selection  of  the  highest  coefficients. 

The  remaining  procedure  is  as  in  Steps  4  through  /  of  Case  1,  except  that 

the  subimages  are  of  size  32  x  64  and  the  feature  vector  V.  .  of  subimage 

$.  .  is  of  dimension  128  as  in  Step  3. 

"1,3 

Case  4:  In  all  the  cases  considered  above,  digital  image  registration  is 
accomplished  by  matching  feature  vectors  whose  elements  frrm  a 
set  of  predetermined  number  of  highest  transform  coefficients. 
Another  criterion  for  the  selection  of  features  for  matching 
purposes  is  to  select  those  transform  coefficients  which  have  the 
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largest  variation  or  variance  ac  jss  image  classes.  Observa¬ 
tion  of  Figures  2-16  through  2-23  indicates  that  most  of  the 
transform  coefficients  with  largest  variance  are  the  low  sequency 
components.  Therefore,  this  case  is  an  effort  to  accomplish 
image  registration  using  features  consisting  of  the  low  sequency 
components. 

Reference  arrays  of  size  32  x  32  are  chosen  from  the  center 
of  the  field  of  view  of  the  reduced  high  resolution  FUR  image. 
The  feature  vector  of  the  reference  image,  V^,  is  formed  by 
selecting  the  63  coefficients  from  the  8x8  array  at  the  top 
left  corner  of  the  transformed  image  excluding  the  (1,1)  element, 
which  represents  the  mean  of  the  pixel  values  of  W.  This  8  x  8 
array  of  transform  coefficients  can  easily  be  computed  without 
actually  computing  the  remaining  coefficients  by  first  reducing 
the  resolution  of  the  reference  array  of  size  32  x  32  by  a  fac¬ 
tor  of  4  to  form  a  reduced  reference  array  of  size  8x8  and  then 
computing  the  Haar  transform  of  the  reduced  array.  The  63-dimen¬ 
sional  feature  vector,  V.  .,  for  each  subimage  S.  .  of  size  32  x 

^  >  J  1,3 

32  from  the  search  area  of  size  120  x  240  (extracted  from  the 
406  x  313  reduced  low  resolution  FLIR  image)  is  also  computed  in 
a  similar  manner.  The  remaining  procedure  of  accomplishing  image 
registration  is  as  in  Case  1. 

Case  5:  Tins  case  is  similar  to  Case  4  except  that  rationalized  Haar 
transform  (RHT)  is  used  to  compute  the  transform  coefficients. 


80 


Case  6:  Here  only  a  subset  of  32  coefficients  from  the  63  coefficients 
selected  in  Case  5,  is  u  >ed  to  form  the  feature  vectors.  The 
rest  of  the  procedure  is  as  in  Case  5.  The  feature  vectors  of 
the  reference  and  s.ubimagcs  of  the  search  area  are  formed  as 
follows. 

1.  An  8  x  8  array  of  transform-coefficients  at  the  upper  left 
corner  of  the  transformed  32  x  32  digital  image  is  computed 
as  explained  in  Case  4. 

2.  Sixteen  coefficients  of  the  4x4  array  at  the  bottom  left 
corner  arid  sixteen  coefficients  of  the  4x4  array  at  the 
upper  right  corner  of  the  8x8  array  of  transform  coeffi¬ 
cients  computed  in  Step  1  are  used  to  form  the  32-dimensional 
feature  vector. 

Case  7:  Reference  arrays  of  size  32  x  64  are  chosen  from  the  center  ol 

the  field  of  view  of  the  reduced  high  resolution  FLIR  image  and 

are  processed  as  two  separate  blocks  of  size  32  x  32  each.  The 
registration  algorithm  is  described  below. 

Step  1:  The  selected  reference  image  of  size  32  x  64  is  divided  into  two 
subimages  and  of  size  32  x  32  each. 

Step  2:  The  resolution  of  W.  and  is  reduced  by  a  factor  of  4  arid  the 

Haar  transform  of  each  reduced  subimage  is  computed  to  form  the 

reduced  transformed  subimages  and  of  sue  8x8  each. 

Step  3:  A  1 26-diinensional  feature  vector,  V  ,  consisting  of  the  coeffi¬ 
cients  of  W-j^  and  exclud  r.g  the  (1,1)  elements,  is  formed. 
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Step  4:  For  each  subimage  S.  .  from  the  search  area  of  size  120  x  240, 

a  126-dimensional  feature  vector  V.  .is  computed  according  to 

i » J 

Steps  1  through  3. 

Step  5:  The  feature  vectors  V..  and  V.  .  are  correlated  to  determine  the 

•  W  i,g 

true  registration  point. 

Case  8:  Similar  to  Case  7  except  that  in  Step  2,  the  modified  Haar  trans¬ 
form  (MHT)  is  used  to  compute  the  transform  coefficients. 

Case  9:  Similar  to  Case  7  except  that  in  Step  2  the  RHT  is  used  to  com¬ 
pute  the  transform  coefficients. 

Case  10:  As  in  Case  7,  with  the  following  changes  in  Steps  3  and  4. 

Step  3:  Sixteen  coefficients  of  the  4x4  array  at  the  bottom  left  corner 
and  sixteen  coefficients  of  the  4x4  array  at  the  top  right 
corner  are  selected  from  both  and  W^R  to  form  a  64-dimen¬ 
sional  feature  vector  V... 

W 

Step  4:  For  each  subimage  S.  .  from  the  search  area  of  size  120  x  240,  a 

1  7  J 

64-dimensional  feature  vector  V.  ■  is  computed  according  to 

'  >J 

Steps  1  through  3. 

Simulation  results  of  the  feature  matching  registration  technique 
using  Haar,  modified  and  rationalized  Haar  coefficients  for  all  the  cases 
considered  above  are  presented  in  Table  12  and  summarized  in  Table  13. 

From  observation  of  the  results  summarized  in  Table  13,  the  following 
conclusions  can  be  made.  The  modified  Haar  transform  has  the  poorest 
performance  and  is  not  suitable  for  extracting  features  for  matching  pur¬ 
poses.  One  possible  reason  for  its  poor  performance  is  the  suppression 
of  the  high  sequency  components  of  .he  transformed  image,  which  represent 
most  of  the  edge  information.  The  RHT,  which  has  an  enhancing  effect  on 
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Table  12.  Continued 


Scenes 


Case  Case  Case 
8  9  10 


Target 


Portion  A 


Target  §2 


Position  B 


Position  C 


Target  #3 


Target  #4 


Target  #1 


Target  #2 


Target  #3 


Target  #4 


Taraet  #1 


Target  #2 


Target  #3 


Target  #4 


Target 


Position  D 


Target  #2 


Target  #3 


Target  #4 


r~ 


HF 


Table  13.  Summary  of  results  in  Table  12. 
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the  edges,  performed  better  than  the  Haar  and  modified  Haar  transforms. 

The  superior  performance  of  RHT  makes  it  more  suitable  for  selection  of 
features  for  matching  purposes  for  the  type  of  scenes  considered.  In  ad¬ 
dition,  since  the  computation  of  RHT  does  not  involve  multiplication  by 
irrational  numbers,  it  can  easily  be  implemented  with  only  adders  and 
subtractors. 

I.  Digital  Image  Registration  Using  Wal sh/Hadamard  Transform 

Walsh  and  Hadamard  transforms,  like  the  Haar  transform,  are  ortho¬ 
gonal  transforms  which  convert  a  digital  image  in  the  "spatial -domain" 
into  a  "transform-domain"  image  containing  the  same  information.  N-di- 
mensional  Walsh  and  Hadamard  transforms,  where  N  is  an  integer  power  of 
two,  possess  a  fast  computational  algorithm  and  differ  only  in  the  order 
in  which  the  rows  and  columns  appear  in  the  transform  matrices.  Since 
the  order  in  which  the  rows  and  columns  appear  is  not  important  for  image 
matching  applications,  the  term  Wal sh/Hadamard  transform  (WHT)  is  used  to 
denote  either  transform.  The  rows  of  the  Wal sh/Hadamard  matrix  are  gen¬ 
erally  ordered  with  increasing  number  of  sign  changes  or  sequency.  The 
sequency  ordered  8x8  Wal sh/Hadamard  matrix  is  shown  in  Figure  2-2 7,  and 
the  two  dimensional  basis  plane  for  an  8  x  8  WHT  is  shown  in  Figure  2-28. 
Wal sh/Hadamard  transform  can  be  implemented  with  only  additions  and  sub¬ 
tractions,  making  it  attractive  for  efficient  hardware  implementation. 

An  effort  is  made  to  study  the  performance  and  accuracy  of  the  fea¬ 
ture  matching  technique  based  on  the  WHT  coefficients.  A  few  simulations 
were  run  for  a  reference  array  of  size  32  x  32  using  a  64-dirnensional 
feature  vector.  In  order  to  reduce  simulation  computation  time,  a  subarray 
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Figure  2-27.  Sequency  ordered  8x8  Wal sh/Hadamard  matrix. 

of  size  130  x  180  which  includes  the  target  area  of  interest  in  extracted 

from  the  416  x  313  reduced  low  resolution  FUR  image  (search  area).  The 

registration  algorithm  is  described  below. 

Algorithm 

Step  1:  A  reference  image  W  of  size  32  x  32  is  selected  from  the  center 
of  the  field  of  view  of  the  reduced  high  resolution  FUR  image. 

Step  2:  The  resolution  of  W  is  reduced  by  a  factor  of  2  to  yield  a  re¬ 
duced  reference  image  of  size  16  x  16,  and  the  two  dimension¬ 
al  WHT  is  computed  by  using  a  fast  Wal sh/Hadamard  transform. 

Step  3:  A  64-dimensional  feature  vector,  V^,  consisting  of  the  64  highest 
coefficients  of  the  transformed  image  W^  is  formed  after  setting 
the  (1,1)  element  of  W^  to  zero.  The  (1,1)  element  represents 
the  mean  of  the  pixel  values  of  W. 

Step  4:  The  WHT  of'  each  subimage  S.  .  of  size  32  x  32  from  the  search 

T  ,  J 

area  is  computed  according  to  Step  2  and  a  64-dimensional  feature 


i  mens 
•ansf' 
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vector  V .  -  is  formed  for  each  subimage  in  such  a  way  that  the 
1  »  3 

elements  of  V.  .  and  V. ,  have  common  locations  in  the  transformed 
i  ,j  W 

images. 

Step  5:  The  feature  vectors  Vw  and  \J.  .  are  correlated  to  find  the  true 
registration  point. 

For  the  purpose  of  comparison,  simulations  were  also  run  by  computing 
Haar  transform  instead  of  WIT  and  the  detailed  results  of  the  simulation 
are  given  in  the  monthly  report  of  August  1981.  A  summary  of  these  re¬ 
sults  is  presented  in  Tables  14  and  15.  The  results  indicate  that  the 
'laar  transform  performs  better  than  the  WHT. 

Another  criterion  for  selection  of  features  for  matching  purposes  is 
to  select  those  transform  coefficients  which  have  the  largest  variation  or 
variance  across  image  classes.  In  order  to  determine  the  coefficients 
with  large  variances,  the  variance  of  each  transform  coefficient  of  the 
Walsh/ Hadamard  transformed  images  of  size  32  x  32  is  computed  for  both  low 
and  high  resolution  FLIR  images  over  the  1600  transformed  images  which  are 
selected  as  outlined  in  a  section  on  the  image  registration  technique 
using  Flaar  transform.  The  variance  distribution  and  the  location  of  the 
largest  variance  transform  coefficients  are  essentially  the  same  as  that 
obtained  using  the  Haar  transform. 

Several  simulations  were  also  run  with  feature  vectors  consisting  of 
the  low  sequency  components  (coefficients  with  large  variance)  of  the 
Walsh/Hadamard  transformed  images.  Even  though  WHT  performed  comparably 
with  the  Haar  transform,  the  performance  and  registration  accuracy  of  the 
rationalized  Haar  transform  were  superior.  It  should  be  noted  that  each 
Walsh  transform  coefficient  is  a  function  of  all  the  pixels  in  the 
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Table  14.  Simulation  results  of  the  feature  matching  technique 
using  WHT  and  Haar  transform  coefficients. 


Scenes  WHT  Haar 


Target  §1  1  1 

Target  #2  X  X 

Position  A 

Target  1!  3  X  X 

Target  #4  X  X 


Target  #1  X  X 

Target  #2  1  1 

Position  B 

Target  #3  2  1 

Target  #4  X  4 


Target  #1  1  1 

Target  #2  X  X 

Position  C 

Target  #3  1  T 

Target  #4  1  1 


Target  #1  X  X 

Target  #2  X  X 

Position  D 

Target  #3  X  X 

Target  #4  X  3 
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Table  15.  Summary  of  results  in  Table  14. 


Method 

Registration  Point 
is  the  First  Peak 

Registration  Point 
is  2nd,  3rd  or 

4th  peak 

Registration  Point 
is  not  in  first 
four  peaks 

WHT 

5 

1 

10 

Haar 

6 

2 

8 
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original  image  (i.e.,  global),  where  as  only  the  first  four  coefficients 
in  the  Haar  case  are  global.  In  other  words,  the  Haar  transform  provides 
a  domain  which  is  sensitive  to  both  global  and  local  variations  of  the 
pixel s  in  an  image. 

The  number  of  additions  and  multiplications  required  to  implement 
the  feature  matching  registration  algorithms  based  on  correlation  of  ad¬ 
jacent  pixels,  difference  between  adjacent  pixels  and  the  Haar  transform 
for  a  search  area  cf  size  240  x  256  and  a  reference  image  (window)  of 
size  32  x  32  are  derived  in  the  previous  sections.  The  total  number  of 
equivalent  real  additions  required  to  implement  correlation,  SSDA,  MHT  or 
RHT,  correlation  of  adjacent  pixels  and  the  new  feature  matching  method 
are  compared  in  Figure  2-29.  Iri  all  the  cases,  the  process  of  deciding 
the  true  registration  point  is  based  on  a  similarity  measure  which  is  ex¬ 
pressed  in  terms  of  the  correlation  of  the  64-dimensional  feature  vectors 
of  the  window  and  the  subimages  of  the  search  area. 


18.0x10 


Figure  2-29.  Comparison  of  the  arithmetic  requirements  of  various  image  registration  algorithms 


III.  HARDWARE  IMPLEMENTATION 


Simulation  results  of  various  registration  algorithms  presented  in 
Chapter  II  indicate  that  the  performance  of  the  moments  method  and  the 
method  based  on  intraset  and  interset  distances  is  unsatisfactory  for  the 
type  of  images  used  for  simulation.  The  other  methods  seem  to  have  com¬ 
parable  performance  to  the  binary  correlation  algorithm.  However,  to  ac¬ 
complish  multiple  image  registration  in  real  time  or  almost  so,  image 
matching  algorithms  must  be  implemented  using  fast  hardware.  An  algorithm 
which  is  computationally  efficient  for  software  implementation  may  be 
complex  or  even  not  feasible  for  hardware  implementation.  Simple  schemat¬ 
ics  shown  in  Figures  3-1  through  3-8  are  used  for  comparison  of  the  com¬ 
plexity  of  hardwarL  implementation  of  the  various  methods.  Since  the 
moments  method  and  the  method  based  on  intraset  and  interset  distances 
are  not  accurate  for  the  kind  of  imagery  considered  in  this  work,  their 
hardware  implementations  are  not  discussed, 

A.  Standard  Cross  Correlation  Algorithm 

No  simulations  are  run  using  the  standard  correlation  algorithm. 

However,  the  standard  correlation  algorithm  is  the  most  popular  and  proven 

method  of  accompli shing  digital  image  registration,  and  therefore  its 

hardware  implementation  is  discussed  in  this  section.  The  normalized 

cross  correlation,  R, (i,j),  between  the  window  W  and  subimage  S.  ■  is 

K  k  T  ,J 

given  by. 
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for 


K  L 

T.  T.  W  (v.,ni)S.  .(fc,m) 
U=1  ni-1  k  1,] 


K  L  o 
l  l  W  (<t,m)l 
£=1  m=l 


'K  L  ? 
l  >.  S  .()t,m) 

3=1  m=l  1>J 


1  <  i  <  M-K+l  and  1  <_  j  <  N-L+l 


(3-1) 


Let 


Ai,j 


K  L 

r.  i  W.  (e.,m)S.  ,U,m) 
a=l  m=l  K  1,J 


K  L  2 

B.  =  z  z  WfU»m) 
K  jl=1  m=l  K 


(3-2) 


(3-3) 


K  L  ? 

C.  i=  £  E  Sj  .(a.m) 

£=1  m=l  11 J 


(3-4) 


Then 


for 


R?d.j)  =  oV— » 

k  BkLi,j 

1  _<  i  _<  M-K+l  and  1  £  j  <  N-L+l 


Ak  , 


(3-5) 


The  schematic  shown  in  Figure  3-1  computes  ( i ,  j )  in  four  stages. 

2  2 

Stage  1:  In  the  first  stage,  S.  .(n.m),  W,  (x. ,m )  and  the  product 

S.  . (x.,m)W.  (t,m)  are  computed  for  n  =  1 ,  2,  3,  K  and  m  = 

1  >  J  K 

1,  2,  ....  L.  This  requires  a  total  of  3KL  two-input  multi¬ 
pliers  (three  multipliers  for  each  pixel  pair). 

Tk  K  L 

Stage  2:  In  the  second  stage,  wA.  .  =  z  z  S.  .(t,m)W,  (n,m) , 

£-1  m=l  1,J  K 

K  L  ?  K  L  2 

B.  s  z  I  W,(£,m)  and  C.  .  =  z  Z  S':  . (£,m)  are  computed 


«.=!  m=l 


llJ  £=1  m=l  1,J 


using  three  KL-input  adders. 


y* 

3 


KL- input 
adder 
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gure  3-1.  Schematic  for  correlation  method. 
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k 

Stage  3:  In  the  third  stage,  two  two-input  multipliers  compute  A.  .  and 

1  *J 

the  product  B,  C .  - . 

K  1  ,  J 

9 

Stage  4:  Finally,  a  divider  computes  R£(i,j)  using  the  results  obtained 
in  Stage  3. 

Therefore,  3KL+2  two-input  multipliers,  three  KL-input  adders  and 

2 

one  divider  are  required  to  compute  R.(i,j)  from  S.  .  and  W,  . 

K  1  >J  K 


B.  Binary  Correlation  Algorithm 

If  the  windows  and  search  area  are  preprocessed  and  transformed  to 

binary  images,  the  correlation  algorithm  can  be  implemented  using  simple 

digital  hardware.  Several  algorithms  for  transforming  digital  images  to 

binary  images  are  given  it.  [5].  When  the  window  and  search  area  are 

in  binary  form,  the  correlation  between  W,  and  S.  .is  defined  as 

K  1  ,J 


}/  ^ 

R  (i,j)  =  KL  -  z  i  |S,  At, m)  -  WU,m}|  ,  (3-6) 

K  *-l  m=l 


for  1  <  i  <  M-K+l  and  1  ^  j  <  N-L+l 


A  schematic  for  the  widely  used  exclusive-OR  gate  implementation  of  the 
binary  correlation  algorithm  given  in  Equation  3-6  is  shown  in  Figure  3-2. 
From  the  schematic,  it  is  clear  that  the  hardware  implementation  requires 
KL  exclusive-OR  gates  and  one  (KL+1)-input  adder.  Detailed  description 
of  this  binary  correlator  can  be  found  in  [2]. 


C.  Correlation  of  Adjacent  Pixels  Method 
An  interesting  property  of  the  image  registration  algorithm  based  on 
correlation  of  adjacent  pixels  is  the  possibility  of  computing  the  fea¬ 
ture  vector  of  the  subimage  S.  ...  or  S.,.  .  from  the  feature  vector  of 

y  i,j+l  i+l,  J 

S.  .  with  very  few  arithmetic  operations.  The  above  property  makes  this 
1 
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method  extremely  suitable  for  real-time  implementation.  A  schematic  for 
a  possible  real-time  hardware  implementation  is  given  in  Figures  3-3 
through  3-5.  It  is  assumed  that  the  windows  and  search  area  have  been 
preprocessed  using  an  algorithm  which  replaces  each  pixel  by  its  value 
minus  the  average  pixel  value  of  the  K1XLI  array  centered  about  it. 

This  compensates  for  the  difference  in  d.c.  bias  between  the  two  imaging 
systems. 

Consider  a  KxN  array  of  shift  registers  as  shown  in  Figure  3-3  where 
K  is  the  number  of  rows  in  the  reference  image  and  N  is  the  number  of 
active  pixels  in  each  row  of  the  preprocessed  search  area.  At  every  sam¬ 
pling  instant  a  new  pixel  of  the  MxN  search  area  enters  the  array  of 
shift  registers  at  the  top  left  corner  as  shown  in  Figure  3-3.  During 
the  first  sampling  period,  S (1 ,1 )  enters  the  shift  register  SR^  ^ .  During 
the  second  sampling  period,  S (1 ,1 )  is  shifted  from  Sf^  ^  to  SR^  2  and 
S(l,2)  enters  SR.  .  and  so  on.  Therefore,  after  N  sampling  periods,  the 

«  5  - 

first  row  of  shift  registers  has  the  N  pixels  of  the  first  row  of  the 
search  area  S.  During  the  N+lth  sampling  period,  S(l,l)  is  shifted  from 
SR.j  ^  to  SR2  the  contents  of  SR^  .  is  shifted  into  SR^  ^  for  1=1, 

2,  ...,  N-l  and  S(2,l)  enters  SR^  The  contents  of  various  shift  reg¬ 
isters  at  the  end  of  KN  sampling  periods  is  shown  in  Figure  3-3.  During 
the  next  sampliny  period  S(l,l)  goes  out  of  the  array  and  $(K+1,1)  enters 
the  array  of  shift  registers  and  so  on.  The  computation  of  feature  vec¬ 
tors  for  all  KxL  subimages  of  the  search  area  as  they  shift  through  the 
array  of  shift  registers  is  presented  next. 
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1 .  Row  Feature  Extractor 

Consider  computation  of  the  normalized  correlation  of  adjacent  pix¬ 
els  in  the  first  row  of  the  subimage  S.  .  using  the  row  feature  extractor 

>  y  I 

shown  in  Figure  3-4.  Multipliers  Ml  and  M2,  accumulators  A1  and  A2,  and 
shi t t  registers  D1  and  02  are  connected  to  SR^  and  SR^  N  as  in  Fig¬ 
ure  3-4.  In  order  to  compute  the  correlation  of  adjacent  pixels  in  row  1 
of  the  subimage  all  shift  registers  and  accumulators  shown  in  Fig¬ 

ure  3-4  are  first  initialized  to  zero.  The  search  area  S  is  gated  into 
the  KxN  array  of  shift  registers  as  described  before.  At  the  end  of 
KN-1  sampling  periods,  the  contents  of  shift  registers  in  row  K  of  the 
array  are  given  in  Figure  3-5.  At  this  point,  M2,  D2,  and  A2  are  enabled. 

p 

It  is  easy  to  see  that  at  the  end  of  KN-1  sampling  periods  S  (1,1),  which 
appears  at  the  output  of  M2,  is  added  to  the  previous  value  in  A2  (zero). 
The  content  of  D2  (L),  which  is  zero  for  this  case,  is  then  subtracted 
from  the  output  of  M2.  Therefore, 

Y  =  S2(l,l)  =  (1,1) 

The  shift  register  D2(l)  now  contains  j(l,l).  Since  Ml,  D1  and  A1  are 

not  enabled  at  the  end  of  KN-1  sampling  periods,  X  is  zero.  At  this 

point.  Ml,  Dl  and  A1  are  enabled.  At  the  end  of  the  next  sampling  period 

(KNth  sampling  period),  SR^  N  and  SR^  ^  ^  contain  S(l,l)  and  S(l,2), 

respectively.  Therefore,  S(l,l)  S(l,2)  appears  at  the  output  of  Ml  and 
2 

S  (1,2)  appears  at  the  output  of  M2.  The  content  of  D2(i)  is  moved  into 
D2(i+1)  for  i  =1,  2,  ...,  L-l ,  and  the  output  cf  M2  is  gated  into  D2(l). 
The  content  of  D2(L)  is  shifted  out  of  D2  and  the  new  value  in  A2  is  the 
old  value  minus  the  contents  of  D2(L)  plus  the  value  at  the  output  of  M2. 
Therefore, 


104 


H 

SO  ,N-1) 

S ( 1  ,N-2) 

•  •  *  j  S(l,l) 

0 

— 

S\,l 

SRK,2 

srk,m-i 

srk,n  _ 

Figure  3-5.  Contents  of  shift  registers  in  row  K  at  the 
end  of  KN-1  sampling  periods. 


Y  =  S2(l,l)  +  S2(l  ,2)  -  0  (3-7) 

a  Sl,l^’^  +  Sl,l^1,2^ 

Similarly, 

X  =  0  +  S(l,l)  S(1 ,2)  -  0  (3-8) 

=  SU1(1,1)  S1  (1,2) 

At  the  end  of  KN+1  sampling  periods, 

Y  =  S2(l ,1 )  +  S2{1 ,2)  +  S2(l ,3)  ,  (3-9) 

and  X  =  S(l,l)  S(l,2)  +  S(1 ,2)  S(1 ,3)  (3-10) 


Also,  D1  (1 )  and  Dl(2)  contain  S(l,2)  S(l,3)  ard  S(l,l)  S(l,2),  respectively. 
Similarly,  02(1),  D2(2),  and  D2(3)  contain  S2(l,3),  S2(l,2)  and  S2 ( 1,1), 
respectively.  During  every  sampling  period  a  product  term  is  added  to  A1 
and  a  square  term  is  added  to  A2.  The  values  in  D1  and  D2  are  shifted 
left  one  place  and  new  values  enter  Dl(l)  and  D2(l).  After  KN+L-2  sampl¬ 
ing  periods. 


L-l  L-l 

X  =  2  S(l,j)  S(l,j+1)  =  z  S,  ,(1,3)  S,  ,0,3+1)  ,  (3-11) 

j=l  3=1  1,1  1,1 

L  2  2 

Y  =  E  S^l.j)  =  Z  S,  ,0,3) 

3=1  3-1 


(3-12) 
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L-l 

2  S,  ,(1,3)  S,  .(1  ,j+l ) 

3  =  1  '  ’ ' _ _ _ 

and  R  =  L  9  (3-13) 

1  I  Sf  ,(1,3) 

3=1  ’ 


From  Equation  (3-13),  it  is  clear  that  R-j  is  the  first  component  of  the 
feature  vector  of  the  subimage  ^ .  At  this  point,  the  shift  registers 

D2(L)  and  Dl(L-l)  contain  S2(l,l)  and  S (1 ,1 )  S(l,2),  respectively.  During 
the  next  sampling  period,  (KN+L-1 )th  sampling  period,  S ( 1 ,1 )  S ( 1 , 2 ) 
shifts  out  of  Dl(L-l)  and  S2(l,1)  shifts  out  of  D2(L).  S(1,L)  S(1,L+1) 

appears  at  the  output  of  Ml  and  S  (1 ,L+1 )  appears  at  the  output  of  M2  and 
therefore, 


and 


L-l 
X  =  z 
3=1 

L-l 
=  z 
j=l 


L 

Y  -  z 
j=l 

L 

=  Z 

j=l 


S(l,j)  S(l,3+1)  +  S(1,L)  S(1  ,L+1 )  -  S(l.l)  S(l,2) 
S1)2(U)  s1>2(1,j+l)  ,  (3-14) 

S2(U)  +  SZ(1,L+1)  -  S2(l,l) 

S2j2(l,j)  ,  (3-15) 

L-l 


z  S1  2(1 ,j)  S1  2(1 ,j+l ) 
—  J  ^  _ 


2  S2  «(1,3) 
j  =  l  1,2 


(3-16) 


From  Equations  (3-16),  it  is  clear  that  R^  is  the  first  component  of  the 
feature  vector  of  the  subimage  ^ •  Similarly,  at  the  end  of  the  next 


I  Ob 

sampling  period,  the  first  component  of  the  feature  vector  of  S-  ,  ap- 

I  ,o 

pears  at  the  output  of  the  divider  and  so  on. 


Z .  Column  Feature  Extractor 

Consider  computation  of  the  normalized  correlation  between  adjacent 
pixels  in  the  first  column  of  the  subimage  S-j  1  using  the  column  feature 
extractor  shown  in  Figure  3-6.  The  SR ^  N  through  SR^  ^  are  the  shift 
registers  shown  in  Figure  3-3.  ,  M2>  ...»  ^ ,  N-j ,  N2>  ....  NK  are 

two-input  multipliers.  At  the  end  of  the  KNth  sampling  period,  the  con¬ 
tents  of  various  shift  registers  of  the  KxN  array  will  be  as  shown  in 
Figure  3-3.  The  first  column  of  ^  will  be  in  the  Nth  column  of  the 
KxN  array  of  shift  registers.  From  Figure  3-6,  it  is  clear  that 


K-l 

X  -  z  S^'i.l)  Sj^O+l.l)  ,  (3-17) 


Y=  Z  Sf  ,(1,1)  ,  (3-18) 

i=l  1,1 


K-l 

z  S,  ,  (i,l )  S.  ,(1+1,1) 

and  ci  =  y  - —  (3-19) 

2  S,2  ,(1,1) 

i=l  1,1 


C1  given  by  Equation  (3-19)  is  the  normalized  correlation  between  the  ad¬ 
jacent  pixels  of  the  first  column  of  ^  and  is  the  (K+l )th  component 

of  the  feature  vector  of  the  subimage  S..  ,.  During  the  next  sampling 

'  >  • 

period,  the  second  column  of  1  will  move  into  the  Nth  column  of  the 
KxN  array  of  shift  registers  shown  in  Figure  3-3  and  therefore  the  new  C^ 

will  be  the  (K+2)th  component  of  the  feature  vector  of  S.  An  efficient 

*  >  1 


nput  adder 
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feature  extractor. 
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and  feasible  way  of  implementing  the  correlation  of  adjacent  pixels  method 
using  K  row  feature  extractors  and  one  column  feature  extractor  is  pre¬ 
sented  next. 


2 .  Real-time  Implementation 

The  schematic  for  the  real-time  hardware  implementation  of  correla¬ 
tion  of  adjacent  pixels  method  using  K  row  feature  extractors  and  a  col¬ 
umn  feature  extractor  is  given  in  Figure  3-7.  From  previous  discussion 
of  the  operation  of  feature  extractors,  it  is  evident  that  at  the  end  of 
the  (KN-l)th  sampling  period, 

R.  =  0  ,  for  1  £  i  <  K  ,  (3-20) 


and  the  correlation  of  adjacent  pixels  in  the  first  column  of  ^  ap¬ 
pearing  at  the  output  of  the  column  feature  extractor  moves  into  SR^.  At 
the  eno  of  KNth  sampling  period. 


R,  = 


Sjj  ( 1 » 1 )  j  ( i ,  2 ) 


1  S^-,  (i  >  I )  +  S^-j  (i  >2) 


,  for  1  <  i  <  K 


(3-21) 


During  the  above  sampling  period,  contents  of  SR..+1  moves  into  SR^  for 
1  £  i  ^  l.-l ,  contents  of  SR j  is  lost  and  the  correlation  of  adjacent  pix¬ 
els  iri  column  2  of  Sj  ^  moves  into  SR^.  Finally,  at  the  end  of  the 
( KN+L-2 ) th  sampling  period, 

L-l 

t  S1  i  (i ,  j )  S^O.j+l) 

R,  =  ^ - 1 - .  for  1  <  i  <_  K  ,  (3-22) 

T  Sf.ti.j) 

3=1  ’’ 


Real-time  implementation  of  correlation  of  adjacent  pixels  method 
and  the  new  feature  matching  method. 


no 


K-l 

z  j  (i  ,j )  i  (i+1 ,  j ) 

and  C.  =  — — —n - ,  for  1  <  j  <  I.  .  (3-23) 

J  k  ?  _  - 

£  s! 

i=l  1,1 

From  Equations  (3-22)  and  (3-23),  it  is  clear  that  [R-j  Ro  •••  C-j 
...  Cl]T  is  the  feature  vector  of  During  the  next,  sampling  period, 

[Rj  R^  •••  Rk  C2  ...  Cl]T  will  be  the  feature  vector  of  S-j  ^  and  s0  cn- 
In  general,  at  the  end  of  the  (KN+L+j-3)th  sampling  period  [R^  R  . .. 

R^  Ci  ...  CLf  will  be  the  feature  vector  of  the  subimage  for 

1  <  j  <  N-L+l .  This  completes  the  computation  of  feature  vectors  for 
subimages  $.  .  for  1  <  j  <  N-L+l.  At  this  poi.it,  all  the  K  row  feature 

•  9  J 

extractors  and  the  column  feature  extractors  are  re-initialized  and  the 

above  procedure  is  repeated  to  compute  the  feature  vectors  of  S~  .  for 

*  »0 

i  <  j  <  N-L+l.  The  above  procedure  is  repeated  until  feature  vectors  of 
all  subimages  are  computed. 

4.  Matching  of  Feature  Vectors 

A  possible  hardware  implementation  for  the  real-time  computation  of 
feature  vectors  of  the  subimages  of  the  search  area  has  been  described  in 
detail.  A  similar  unit  consisting  of  a  KxL  array  of  shift  registers,  K 
row  feature  extractors  and  a  column  feature  extractor,  can  be  used  to  com¬ 
pute  the  feature  vectors  of  the  n  windows  sequentially  or  n  such  units 
can  be  used  in  parallel  to  compute  the  n  feature  vectors  simultaneously. 

In  order  to  accomplish  multiple  image  registration,  the  feature  vector  of 
each  window  must  be  matched  with  the  feature  vectors  of  the  subimages  of 

5.  One  way  of  accomplishing  the  above  is  shown  in  Figure  3-8  and  is  self 
explanatory. 


ure  3-8.  Correlation  of  feature  vectors. 
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D.  Hardware  Requirements  of  the  Standard  Correlation  and 
Correlation  of  Adjacent  Pixels  Method 

A  schematic  for  real-time  implementation  of  the  standard  correlation 
algorithm  is  given  in  Figure  3-1.  For  simultaneous  registration  of  n 
windows,  n  such  units  must  be  used  in  parallel.  C.  •  in  Equation  3-4  is 
independent  of  the  windows  and  need  be  computed  only  once.  iherefore, 

3KL+2  two-input  multipliers,  one  divider  and  three  KL-input  adders  are 
required  to  implement  the  first  correlator.  Each  additional  correlator 
can  be  implemented  using  2KL+2  two-input  multipliers,  one  divider  and  two 
KL-input  adders.  In  other  words,  simultaneous  registration  of  n  windows 
using  the  standard  cross  correlation  method  requires  (2KL+2)n  +  Kl.  two- 
input  multipliers,  (2n+l)  KL-input  adders  and  n  dividers. 

Each  of  the  K  row  feature  extractors  shown  in  Figure  3-7  requires  2 
multipliers,  2  accumulators,  one  divider  and  2L+1  shift  registers.  The 
column  feature  extractor  requires  2K-1  multipliers,  one  (K-l)-input  adder, 
one  K  input  adder,  one  divider  and  L  shift  registers.  Therefore,  the 
schematic  in  Figure  3-7  requires  a  total  of  4K-1  multipliers,  2K  accumula¬ 
tors,  2KL+K+L  shift  registers  (excluding  the  KxN  array  of  shift  registers), 
one  K-input  adder,  one  (K-l)-input  adder  and  K+l  dividers.  The  feature 
vectors  of  all  the  n  windows  are  computed  first  sequentially,  using  the 
hardware  shown  in  Figure  3-7,  and  then  the  feature  vectors  of  the  subimages 
of  the  search  area  are  computed.  The  matching  of  feature  vectors  shown  in 
Figure  3-8  requires,  {(2K+2l+l)n  +  K+L}  multipliers,  n  dividers  n  square 
root  operators,  and  (2n+l)  (K+L)-input  adders.  The  above  results  are  sum¬ 
marized  in  Table  16.  If  the  search  area  and  windows  are  coarsely  quanti¬ 
zed  to  two  or  three  levels,  multipliers  can  be  implemented  using  exclusive- 
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Table  16.  Hardware  requirements  of  the  standard  correlation 
and  correlation  of  adjacent  pixels  methods. 


Standard  Correlation 
Method 

Correlation  of  Adjacent 
Pixels  Method 

Multipl iers 

(2KL+2)n  +  KL 

{ 4  K-l )  +  ( 2K+2L+1 )n  +  K  +  L 

Dividers 

n 

K  +  1  +  n 

Adders 

2n+l  KL-input 

one  K-input 
one  (K-l ) -input 

2n+l  (K+L)-input 

Shift  Registers 

KM  +  nKL 

KN  +  2KL  +  K  +  l. 

Square  Root 
Operator 

- 

n 

Accumulators 

- 

2K 
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OR  gates  and  the  implementations  of  the  standard  correlation  algorithm 
and  correlation  of  adjacent  pixel  methods  become  much  simpler  to  realize. 

E.  New  Feature  Matching  Method 

In  this  method,  digital  image  registration  is  accomplished  by  using 

features  based  on  adjacent  pixel  differences  along  the  rows  and  columns 

of  a  digital  image.  Once  the  feature  vector  of  the  subimage  S.  .is  com- 

^  >  J 

puted,  the  feature  vector  of  the  next  subimage,  .+^ ,  can  be  computed 
with  very  few  additional  computations.  This  recursive  nature  of  the  com¬ 
putation  of  feature  vectors  makes  the  new  feature  matching  method  effi¬ 
cient  for  real-time  ha> dware  implementation.  A  KxN  array  of  shift  regis¬ 
ters  as  shown  in  Figure  3-3  initially  containing  the  first  K  rows  of  the 
search  area  of  size  MxN  is  used  for  the  purpose  of  extracting  row  and 
column  features.  At  each  sampling  instant,  a  new  pixel  from  the  search 
area  enters  the  shift  register,  SR^  ^ ,  pushing  the  previous  content  of 
SR^  ^  and  the  contents  of  other  registers  by  one  place  to  the  right.  A 
possible  hardware  implementation  of  the  row  and  column  feature  extractors 
is  discussed  first  and  a  schematic  for  real-time  implementation  of  the 
feature  matching  algorithm  is  presented  next. 

1 .  Row  Feature  Extractor 

A  hardware  implementation  of  the  row  feature  extractor  is  shown  in 
Figure  3-9.  The  row  feature  extractor,  connected  to  the  shift  registers 
SR^  ^  and  SR^  ^  ^  in  a  KxN  array  of  shift  registers  which  stores  a  part 
of  t  search  image,  is  used  to  compute  the  row  feature  element  of  the 
first  row  of  the  subimage  of  size  Kxl.  Initially  the  shift  registers  SR^, 
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and  SR2  and  the  accumulators  A.|  and  are  cleared.  The  operation  of  the 
row  feature  extractor  is  similar  in  some  respects  to  the  one  described 
under  correlation  of  adjacent  pixels  method  and  is  briefly  described  be¬ 
low. 

During  the  next  sampling  instant,  the  contents  of  shift  registers, 

SR^  ^  and  SR^  ^  -j,  are  moved  into  the  subtractor  and  the  content  of 
SR^  which  is  S(!,l),  is  loaded  into  the  shift  register  SR,,(L)  and  ac¬ 
cumulator  At  A,,,  the  previous  content  of  A^  (zero)  is  added  to  S ( 1 ,1 ) 
and  the  content  of  SR^O),  which  is  also  zero,  is  subtracted  from  S ( 1 ,1 ) . 
At  every  sampling  period  the  contents  of  the  shift  registers  in  the  KxN 
array  are  shifted  to  the  right  by  one  place,  which  results  in  contents  of 
SR  .  occupying  SR..  ...  During  the  next  sampling  period,  S (1 ,2)  and 

IN*  11-  I  t\  3  |V 

$(1,3)  are  moved  into  the  subtractor  and  at  the  same  time  the  output  of 
the  subtractor  which  is  the  absolute  value  of  the  difference  between 
S(l,l)  and  S(l,2)  is  loaded  into  SR  (L-l )  and  accumulator  A^ .  At  A^ ,  the 
output  of  the  subtractor  is  added  to  the  previous  contents  of  A-j  and  the 
content  of  SR^(l)  is  subtracted  both  of  which  are  zero.  Also  at  this 
sampling  instant,  S ( 1 , 2 )  is  moved  into  SR^L)  and  where  the  previous 
content  of  A^,  which  is  S(l,l),  is  added  to  S(l,2)  and  the  content  of 
SRZ(1)  (zero)  is  subtracted.  After  the  Lth  sampling  instant,  the  quo¬ 
tient  of  the  contents  of  A^  and  A^  will  give  the  row-feature  element 
R^(l,l)  of  the  first  subimage  with  R^  appearing  at  the  output  of 

the  divider  after  the  (l+l)th  sampling  period.  The  feature  element 
R|(l,l)  can  be  expressed  in  terms  of  the  elements  of  the  first  row  of 
the  subimage  ^  according  to  the  equation 


MU) 

Ri(ul)  =  dJtuT 
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(3-24) 


where 


^  (1,1 ) 

L-l 

=  I 

j=l 

|S1 J (l,j)  -  s] J (1 »j+l )  1  r  and 

0^1,1) 

L 

=  >: 

|slsl(i,j)| 

During  the  (L+l)th  sampling  period,  the  content  of  SR^(l),  which  is 

|S(1,1)  -  S(l,2)|,  is  subtracted  from  A^  and  |S(1,L)  -  S(1,L+1)|  is  added 

At  A^,  S(1,L+1)  is  added  to  and  S(l,l)  is  subtracted  from  the  content  of 

A,,,  so  that  the  output  of  the  divider  at  the  (L+2)th  sampling  instant  is 

R| (1 ,2) ,  the  row-feature  element  of  the  subimage  S1  according  to  the 
eolation 


R-|  (1  >2) 


Nt(1  ,1  )-|S1 Q  ,1  )-S]  J  (1 ,2)|  +  |$1 11  (1  .Q-Sj  (1  >L+1 ) 


D1(1,'I)-|S11(1,1)|  +  |S1  sl(l,  L+l )  l 


(3-25) 


Similarly  after  (L+3)  sampling  periods,  Rj(l,3),  the  row  feature  element 
of  the  subimage  -  appears  at  the  output  of  the  divider  and  so  on. 


2 •  Column  Feature  Extractor 

A  possible  hardware  implementation  of  the  column  feature  extractor 
is  shown  ir.  Figure  3-10.  It  is  again  assumed  that  a  part  of  the  search 
area  is  stored  in  the  KxN  array  of  shift  registers  shown  in  Figure  3-3, 
with  the  first  column  of  the  subimage  ^  now  residing  in  the  last  col¬ 
umn  of  the  array  to  which  the  feature  extractor  is  connected.  The  adders 
A,  through  A.,  ,  compute  tiie  adjacent  pixel  differences  of  the  elements  of 
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the  first  column  whose  absolute  values  are  added  together  by  a  ( K-l ) - 
input  adder.  A  K-input  adder  shown  in  Figure  3-10  finds  the  sun;  of  all 
the  pixels  in  the  first  column  of  the  subimage.  The  outputs  of  the  K- 
input  and  (K-l)-input  adders  are  connected  to  a  divider  whose  output,  C-j , 
is  the  column-feature  element  of  the  first  column  of  the  subimage  ^  as 
given  by  the  equation, 

K-l 

Z  |S(i,l)  -  S(i+1 ,1 )  | 

C]  =  — -  (3-26) 

z  |  S  ( i ,  1 )  l 
i=l 

During  the  next  sampling  period,  the  second  column  of  the  subimage  ^ 
will  move  into  the  last  column  of  the  KxN  array  of  shift  registers  and 
the  column  feature  element,  C of  the  second  column  of  the  subimage  will 
be  computed  by  the  column  feature  extractor.  The  process  continues  at 
each  succeeding  sampling  period. 

3 .  Real-time  Hardware  Implementation  of  the 
New  Feature  Hatching  Method 

An  efficient  way  of  implementing  the  new  feature  matching  algorithm 
using  the  row  and  column  feature  extractors  described  in  the  previous 
sections  is  shown  in  Figure  3-7.  The  K  row-feature  extractors  compute 
the  feature  elements  of  the  K  rows  of  the  subimage  in  L  sampling  periods. 
The  column  feature  extractor  computes  the  column-feature  element  of  any 
one  column  of  the  subimage  in  one  sampling  period,  thus  computing  all  the 
l  column-feature  elements  at  the  end  of  L  sampling  periods.  Initially, 
it  takes  L  sampling  periods  to  compute  the  (K+L)-dimensional  feature  vec¬ 
tor  i  of  the  subimage  ^ .  Next  it  takes  only  one  additional  sampling 
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period  to  compute  the  feature  vector  ^  the  suhimage  ^  and  so  on, 
till  all  the  feature  vectors  along  a  row  of  the  search  image  are  computed. 
The  process  is  repeated  to  compute  the  feature  vectors  along  the  other 
rows  of  the  search  area,  reinitializing  (clearing  the  registers  and  accu¬ 
mulators)  the  row-and  column-feature  extractors  every  time  the  computa¬ 
tion  starts  at  the  beginning  of  a  new  row. 

4 •  Multiple  Image  Registration 

Multiple  image  registration  is  accomplished  by  making  a  decision 
based  on  a  similarity  measure  which  is  expressed  in  terms  of  the  correla¬ 
tion  of  the  feature  vectors  of  the  subimages  with  the  feature  vector  of 
each  window.  A  hardware  implementation  used  for  this  purpose  is  illus¬ 
trated  in  Figure  3-8. 

5 .  Hardware  Requirements  of  the  New 
Feature  Matching  Method 

The  complexity  of  the  hardware  required  for  simultaneous  registration 
of  n  windows  using  the  real-time  hardware  shown  in  Figure  3-7  is  given 
below. 

1.  A  row-feature  extractor  requires  one  subtracter  (adder),  one 
divider,  two  accumulators  and  (21+1 )  shift  registers. 

2.  A  column-feature  extractor  requires  (K-l)  adders,  one  divider, 
one  K-input  adder,  one  (K-l)-input  adder  and  L  shift  registers. 

3.  The  correlation  of  feature  vectors  of  subimages  and  windows 
for  matching  purpose  (Figure  3-8)  requires  (2K.+2L+1  )n+K+L 
multipliers,  n  dividers,  n  square  root  operators,  and  (2n+l) 
(K+L)-input  adders. 
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The  total  hardware  requirement  for  the  real-time  implementation  of 
the  feature  matching  multiple  registration  algorithm  as  illustrated  in 
Figure  3-7,  with  K  row-feature  extractors,  one  column-feature  extractor 
and  a  correlator  for  matching  feature  vectors  is  summarized  in  Table  17. 

F.  Haar  Transform  Method 

The  Flaar  transform  can  be  computed  efficiently  using  a  fast  Haar 
transform  (FHT).  A  flow  diagram  for  computing  an  8-point  FHT  is  shown  in 
Figure  2-13.  A  two-dimensional  Haar  transform  is  then  computed  by  re¬ 
peated  application  of  the  one-dimensional  FHT.  A  2-dimensional  transform 
of  a  digital  image  of  size  KxK  is  obtained  by  first  computing  the  Haar 
transform  of  each  of  the  K  rows  (horizontal  transformation)  and  then  com¬ 
puting  the  transform  of  each  of  the  K  columns  (vertical  transformation) 
of  the  resulting  KxK  image  after  horizontal  transformation.  Some  varia¬ 
tions  of  the  basic  Haar  transform,  such  as  the  modified  Haar  transform 
(MHT)  and  the  rationalized  Haar  transform  (RHT)  do  not  involve  multiplica¬ 
tion  by  irrational  numbers  (  /Z)  and  can  be  implemented  with  adders  and 
subtractors  only,  making  them  more  attractive  for  real-time  applications. 
It  has  also  been  shown  in  the  previous  chapter  that  RHT  performs  better 
than  the  Haar  and  the  modified  Haar  transform  in  registering  two  images 
taken  from  different  sensors  for  the  scenes  considered.  A  possible  hard¬ 
ware  implementation  of  the  two-dimensional  Haar  transform  using  both 
pipelined  and  parallel  one-dimensional  Haar  transform  processing  units  is 
discussed  in  the  following  sections. 
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Table  17.  Hardware  requirements  of  the  new 
feature  matching  method. 

Adders  2K-1  2-input  adders 

(subtractors)  2n+l  (K+L)-input  adders 

one  K-input  adder 

one  (K-l )-input  adder 


Multipliers  (2K+2L+l)n  +  K  +  L 


Dividers  K  +  1  +  n 


Shift  registers  (2L+1)K  +  L  +  KN 


Accumulators  2K  +  L 


Square  root  n 

opera  tors 
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1 .  Pipelined  Implementation  of  the  Haar  Transform 

Pipelined  implementation  of  the  8-point  FHT  shown  in  Figure  3-11  is 
best  understood  by  referring  to  the  flow  diagram  shown  in  Figure  2-13. 

The  pipeline-transform  unit  shown  in  the  figure  is  connected  to  the  last 
row  of  the  KxN  array  of  shift  registers  storing  a  part  (K  rows)  of  the 
search  area.  Initially,  the  first  K  rows  of  the  search  area  are  stored 
in  the  array  and  the  Haar  transform  of  the  first  row  of  the  subimage  ^ 
of  size  KxK  for  K  =  8  is  computed  in  the  pipeline-transform  unit  as  fol¬ 
lows.  The  pixels  of  the  first  row  of  ^  enter  the  transform  unit  se¬ 
quentially  one  at  a  time  at  each  sampling  instant.  The  first  pixel  is 
stored  temporarily  in  latch  and  is  presented  to  the  first  input  of 
the  adder/ subtractor  AS^  of  stage  1  during  the  second  sampling  instant 
at  which  time  the  second  pixel  appears  at  the  input  of  the  transform  unit 
and  also  at  the  input  I2  of  AS^ .  At  this  point,  the  sum  and  difference 
of  the  first  two  pixels  are  computed,  storing  the  sum  temporarily  in  the 
shift  register  D^,  while  the  difference  is  outputted  as  one  of  the  trans¬ 
form  coefficients.  The  processing  is  identical  at  each  stage;  the  sum 
and  difference  of  successive  terms  from  the  preceeding  stage  are  computed; 
the  sum  is  propagated  to  the  following  stage,  and  the  difference  is  out¬ 
putted  as  a  transform  coefficient.  It  is  to  be  noted  that  two  shift  reg¬ 
isters  and  are  needed  for  the  second  stage  and  four  shift  registers 
through  D j  are  needed  for  the  third  stage  in  order  to  provide  proper 
inputs  to  the  Adder/Subtractors  for  the  computation  of  transform  coeffi¬ 
cients.  Initially,  a  delay  of  seven  sampling  periods  (neglecting  the  de¬ 
lay  involved  in  the  computation  of  sum  and  differences  in  the  adder/ 
subtractors)  is  needed  to  get  all  the  8  transform  coefficients  of  the 
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first  row  of  the  subimage  S,  During  the  next  sampling  period,  all  8 

I  »  J 

transform  coefficients  of  the  first  row  of  the  subimage  ^  are  available 
at  the  output  since  the  input  sequence  of  pixels  of  the  first  row  of  ^ 
are  the  input  sequence  of  pixels  of  the  first  row  of  ^  shifted  one 
place  to  the  right,  and  differing  only  in  the  last  pixel.  During  the 
next  sampling  instant,  the  coefficients  of  the  first  row  of  the  subimage 
S.j  ^  is  available  at  the  output  and  so  on. 

2.  Parallel  Implementation  of  the  Haar  Transform 

Parallel  implementation  of  the  8-point  fast  Haar  transform  using 
adder/ subtractors  shown  in  Figure  3-12  is  the  hardware  realization  of  the 
FHT  flow  diagram  shown  in  Figure  2-13.  The  operation  of  the  parallel 
transform  unit  is  easily  understood  from  the  flow  diagram  of  the  Haar 
transform.  A  sequence  of  eight  pixels  whose  transform  is  to  be  computed 
is  input  to  the  transform  unit  at  the  same  instant  of  time  and  the  trans¬ 
form-coefficients  are  available  at  the  output  after  a  small  amount  of 
delay  occurring  in  the  adder/subtractors.  This  delay  should  not  exceed 
one  sampling  instant  for  efficient,  real-time  implementation  of  the  fea¬ 
ture  matching  method  using  Haar  coefficients  as  discussed  in  the  next  sec¬ 
tion.  An  increase  in  speed  is  achieved  by  the  parallel  transform  hardware 
at  the  expense  of  an  additional  four  adder/subtractors  over  that  required 
for  the  pipeline  transform  unit  shown  in  Figure  3-11. 

3 .  Real-Time  Hardware  Implementation  of  the 
Haar  Transform  Method 

A  hardware  implementation  of  the  feature  matching  algorithm  based  on 
Haar  coefficients  using  pipelined  and  parallel  Haar  transform  units  des¬ 
cribed  in  the  previous  sections  is  illustrated  in  Figure  3-13.  The  K 
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Real-time  implementation  of  the  feature  matching  method 
based  on  Haar  coefficients. 
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pipelined  Haar  transform  units  (HT  )  connected  to  the  last  column  of 
shift  registers  1  f  the  KxN  array  of  shift  registers,  compute  the  row 
transform  of  the  K  rows  of  the  subimages  of  size  KxK.  Row  transformation 
of  the  first  subimage  requires  seven  sampling  instants,  with  subsequent 
subimages  requiring  only  one  sampling  instant.  The  row  transformed  sub¬ 
image  is  stored  in  a  KxK  array  of  shift  registers.  The  K  parallel  Haar 

transform  units  (HT  )  connected  to  each  column  of  the  KxK  array  of  shift 
P 

registers  storing  the  row  transformed  image,  compute  the  column  transform 
of  the  image  in  one  sampling  instant,  thus  computing  the  two-dimensional 
transform  of  the  original  subimage.  The  KxK  array  of  shift  registers  is 
updated  with  a  >■  ew  row  transformed  subimage  at  every  sampling  instant, 
and  the  paralle'  transform  units  compute  the  column  transform  so  that  the 
two  dimensional  ransform  of  the  new  subimage  is  available  during  the  next 
sampling  instant  Each  time  the  computation  starts  at  the  beginning  of  a 
new  row  of  the  ■.  arch  image,  there  is  a  delay  of  eight  sampling  instants 
before  the  two  d  mensional  transform  of  the  subimage  is  computed. 

4.  Multiple  Image  Registration 

All  or  a  subset  of  the  Haar  transform  coefficients  computed  are  used 
to  form  a  feature  vector  of  a  subimage  or  a  window.  Multiple  image 
istration  is  then  accomplished  by  making  a  decision  based  on  a  similarity 
measure  which  is  expressed  in  terms  of  the  correlation  of  feature  vectors 
of  the  subimages  with  the  feature  vector  of  each  window.  A  hardware  im¬ 
plementation  used  for  this  purpose  is  shown  in  Figure  3-8. 
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5 .  Hardware  Requirements  of  the  Ha ar  Transform  Method 

The  amount  of  hardware  required  for  simultaneous  registration  oi  ,, 
windows  using  the  real-time  hardware  shown  in  Figure  3-13  is  as  follows. 

1.  A  pipelined  Haar- transform  unit  computing  K-point  FHT, 

p 

where  K  is  an  integer  power  of  2  (K  =  2  ),  requires  P 
adder/subtractors  and  (K-l)  shift  registers,  where 
P  =  log2K. 

2.  A  parallel  Haar-transform  unit  requires  (K-l)  adder/ 
subtractors. 

2 

3.  A  KxK  array  (i.e.,  K  )  of  shift  registers  is  needed  to 
store  the  row  transformed  image. 

4.  The  correlation  of  R-dimensional  feature  vectors  of  sub¬ 
images  and  windows  for  matching  purposes  (Figure  3-8) 
requires  (2R+1  )n+R  multipl  iers,  ri  dividers,  n  square 
root  operators,  and  (2n+l )  R-input  adders. 

The  dimension  R  depends  on  the  number  of  Haar  coefficients  selected 
for  matching.  The  total  hardware  requirement  for  real-time  implementation 
of  the  feature  matching  method  using  Haar  coefficients  with  K  pipelined 
transform  units  and  with  K  parallel  transform  units  is  summarized  in  Table 
18. 
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Table  18.  Hardware  requirements  of  the  feature  matching 
method  using  Haar  coefficients. 


Adder/ Subtractors 
Shift  registers 
Mul tipi iers 
Dividers 

Square  root  operators 


K[( K-l )  +  log2K] 
K(K-1 )  +  K2  +  KN 
(2R+1 )n  +  R 


R-input  adders 


(2n+l) 


IV.  HARDWARE  vs.  TIME  TRADEOFFS  FOR  MULTIPLE 
TARGET  HANDOFF  PROBLEM 

A  number  of  feature  matching  algorithms,  their  image  matching  accu¬ 
racies  obtained  through  simulation  and  suggested  hardware  implementations, 
are  given  in  Chapters  II  and  III.  In  this  chapter,  several  methods  of 
accomplishing  multiple  target  handoff  are  presented.  The  trade-off  be¬ 
tween  amount  of  hardware  required  and  total  estimated  time  for  handing 
off  n  targets  to  n  missile  seekers  is  discussed. 

A.  Method  One 

Figure  4-1  with  correlators  2  through  8  not  present  is  a  block  dia¬ 
gram  of  a  single  target  hand-off  correlator  modified  to  hand  off  n  tar¬ 
gets  sequentially.  Initially  the  helicopter  unmasks  and  targets  are 
recognized  and  cued  in  the  acquisition  sensor  field  of  view  either  manu¬ 
ally  or  using  some  automatic  cueing  device  such  as  moving  target  indica¬ 
tor  (MTI)  radar.  MICOM  has  demonstrated  in  its  hand-off  technology  pro¬ 
gram  that  manual  cueing  of  targets  is  too  time  consuming,  exposing  the 
helicopter  to  hostile  fire,  thereby  reducing  the  survivability  of  the 
aircraft.  Target  X  and  Y  positions  within  the  acquisition  sensor  FOV  are 
then  loaded  into  registers  in  the  target  assignment  circuitry. 

As  pointed  out  in  Chapter  1,  the  FOV  of  the  acquisition  sensor  is 
assumed  to  be  larger  than  the  FOV  of  the  missile  seeker  sensors.  Also, 
it  is  assumed  that  the  LOS  of  each  missile  sensor  has  been  prealigned 
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with  the  acquisition  sensor  LOS  such  that  the  field  of  view  of  each 
seeker  is  within  the  FOV  of  the  acquisition  sensor.  The  next 
task  in  the  multiple  target  handoff  problem  is  to  locate  the  missile 
sensor  LOS  positions  within  the  acquisition  sensor  FOV.  In  method  one, 
this  is  done  sequentially  by  loading  a  reference  array  from  the  center 
of  the  FOV  of  the  first  missile  seeker  into  correlator  ft  1  shown  in  Figure 
4-1.  The  adaptive  digitizer,  prefiltering,  edge  detection,  and  slicing 
blocks  in  Figure  4-1  serve  tne  same  functions  as  in  a  one-to-one  handoff 
correlator.  Acquisition  sensor  video  is  then  preprocessed  and  used  as 
the  live  array  (search  array)  by  the  correlator  circuitry.  The  peak  in 
the  correlation  surface  is  the  sensor  X  &  V  LOS  location  within  the  search 
area  and  is  stored  in  a  register  within  the  correlator  unit.  Four  addi¬ 
tional  search  area  fields  are  correlated  with  the  reference  array  and  the 
X  &  Y  I  .OS  locations  are  loaded  into  additional  registers  in  the  correla¬ 
tor.  An  algorithm  is  then  used  to  determine  the  most  probable  LOS  loca¬ 
tion  or  to  determine  if  in  fact  a  valid  registration  has  been  achieved. 
MICOM  has  demonstrated  in  its  handoff  technology  program  that  correlating 
on  five  successive  fields  is  sufficient  to  make  the  valid  match  point 
decision  [11]. 

The  valid  match  point  X  &  Y  values  are  then  used  by  the  target  as¬ 
signment  circuitry  to  assign  a  specific  cued  target  within  the  acquisition 
FOV  to  the  missile  seeker  using  a  nearest  target  to  sensor  LOS  rule.  Er¬ 
ror  signals  aX,  and  aY,  are  then  generated  and  used  to  cause  the  missile 
LOS  to  slew  in  the  direction  of  its  assigned  target.  It  is  assumed  t hat 
slewing  can  be  accomplished  in  0.3  seconds  or  less.  This  time  is 
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obviously  dependent  on  the  line  of  sight  error  magnitude  and  on  the  slew 
rate  of  the  seeker. 

Concurrent  with  the  slewing  of  missile  number  one  LOS  the  identical 
procedure  described  in  the  above  two  paragraphs  is  repeated  using  video 
from  missile  number  two  sensor.  The  process  is  repeated  sequentially 
for  up  to  eight  missiles.  It  is  felt  that  after  the  first  slewing  a 
second  look  at  each  missile  LOS  will  be  necessary  to  verify  that  it  is  on 

the  right  target  and  that  its  LOS  is  close  enough  to  the  target  for  the 

seeker  tracker  to  lock  on.  If  it  is  not,  then  a  second  slewing  will  be 
required. 

In  computing  time  required  for  firing  a  total  of  eight  missiles,  two 
assumptions  have  to  be  made.  First,  in  order  for  the  first  missile  fired 

not  to  blind  the  last  missile  fired,  it  is  assumed  that  the  total  elapsed 

time  between  the  first  and  last  missile  firing  will  not  exceed  1.75  sec¬ 
onds.  Second,  in  order  to  allow  launch  transients  to  die  out  on  a  given 
missile  launcher,  firing  of  two  missiles  from  the  same  launcher  shall  not 
be  closer  than  0.35  seconds.  If  two  launchers  are  available  and  the  mis¬ 
siles  are  fired  in  tandem  from  the  two  launchers,  then  a  missile  can  be 
fired  every  0.175  seconds.  It  should  be  pointed  out  that  these  are  as¬ 
sumptions  with  the  exact  firing  time  constraints  established  in  a  tech¬ 
nology  program. 

Loading  of  each  missile  seeker  video  will  require  2  fields  (one  to 
allow  for  proper  sync  before  loading  of  field  of  data).  Loading  of  the 
first  acquisition  sensor  field  will  also  require  2  fields  with  each  four 
succeeding  fields  loaded  in  real-time.  Assignment  of  targets  is  accom¬ 
plished  by  a  dedicated  processor  and  done  concurrently  with  correlation 


of  the  following  missile  video.  Therefore,  in  the  computation  of  total 
firing  time,  the  assignment  of  targets  is  assumed  to  require  no  addition¬ 
al  time.  Slewing  of  each  missile  LOS  is  assumed  to  be  0.3  seconds  or 
less  and  is  done  concurrently  with  correlation  of  other  missile  seekers 
therefore  requiring  no  added  time.  The  check  of  each  missile  LOS  to  in¬ 
sure  proper  slewing  will  require  eight  fields  (2+2+4).  Any  additional 
slewing  is  assumed  to  be  done  concurrently  with  the  missile  firings.  It 
is  assumed  that  at  least  one  missile  will  not  require  any  additional  LOS 
slewing  and  can  be  fired  immediately  upon  completion  of  the  check. 

From  the  above,  the  first  missile  can  be  fired  in  [8(2+2+4)  +  8] 
(1/60)  =  1.2  seconds.  All  eight  missiles  cart  be  fired  in  1.2  +  7(0.175) 
=  2.43  seconds  if  all  missiles  are  fired  with  a  minimum  time  interval. 

The  advantage  of  this  method  is  that  it  requires  a  minimum  amount 
of  hardware  change  from  a  single  target  handoff  correlator.  The  disad¬ 
vantage  is  that  it  requires  more  time  than  other  methods  because  of  the 
sequential  processing. 

Using  this  method,  the  total  time  required  to  fire  six  missiles  is 
[6(2+2+4)  +  8](l/60)  +  5  (0.175)  =  1.81  seconds  and  the  time  required  to 
fire  four  missiles  is  [4(2+2+4)  +  8](l/60)+  3  (0.175)  =  1.19  seconds. 

The  firing  times  for  the  first  missile  fired,  the  firing  of  eight,  six 
and  four  missiles  are  given  in  Table  4-1  at  the  end  of  this  chapter  for 
comparison  purposes. 

B.  Method  Two 

The  block  diagram  for  method  two  is  shown  in  Figure  4-1  with  all 
eight  correlators  present.  Targets  are  cued  in  the  acquisition  sensor 
F0V  the  same  way  they  were  in  method  one.  A  reference  array  from  the 
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center  of  the  first  seeker  video  FOV  is  loaded  into  correlator  ff\ .  Re¬ 
ference  arrays  from  the  remaining  seven  sensors  are  then  loaded  into 
correlators  two  through  eight  sequentially.  This  will  require  a  maximum 
of  16  video  fields  (two  fields  for  each  reference  to  allow  for  sync  er¬ 
rors).  The  acquisition  sensor  video  is  then  preprocessed  and  used  as 
the  line  (or  search)  array  in  all  eight  correlators  in  parallel.  Again, 
correlation  is  performed  over  five  fields  in  order  to  insure  registra¬ 
tion  accuracy  and  reliability.  A  total  of  six  fields  is  required  for 
this  step.  In  this  step,  LOS  position  results  from  the  five  correlation 
surfaces  are  averaged  after  obvious  non-correlation  fields  are  thrown 
out.  The  X  and  Y  LOS  correlation  results  from  all  eight  sensors  are 
sent  to  the  target  assignment  circuitry  in  parallel  and  targets  are  as¬ 
signed  to  each  missile  seeker  using  a  minimum  total  slew  distance  rule. 
This  step  is  assumed  to  require  1/60  second.  Error  signals  are  then 
generated  thereby  causing  each  missile  LOS  to  slew.  After  waiting  0.3 
seconds  to  allow  for  missile  LOS  slewing,  the  above  steps  are  repeated. 

At  this  point,  it  is  assumed  that  several  missile  LOS  are  aligned  such 
that  their  trackers  can  lock  on  to  their  target  arid  firing  can  occur. 

Any  additional  slewing  of  missiles  which  are  not  aligned  to  the  desired 
accuracy  is  assumed  to  be  done  while  the  first  missiles  are  being  fired 
and  therefore  does  not  require  additional  time. 

From  the  above,  the  first  missile  can  be  fired  in  2(16+6+1 ) (1/60)+ 
0.3  =  1.07  seconds  and  all  eight  missiles  can  be  fired  in  1.07  +  7(0.175) 
=  2.2S  seconds.  The  advantage  of  this  method  is  that  it  offers  a  slight 
total  firing  time  improvement.  The  disadvantage  is  that  eight  full  cor¬ 
relators  are  required. 
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For  comparison  with  method  one,  six  missiles  can  be  fired  in 
2(1 2+6+1 ) (1/60)  +  0.3  +  5  (0.175)  =  1 .81  seconds  and  four  missiles  can 
be  fired  in  2 (8+6+1 ) (1/60)  +  0.3  +  3  (0.175)  =  1.33  seconds. 

C.  Method  Three 

The  block  diagram  for  method  three  is  shown  in  Figure  4-1  with  only 
correlators  1  through  4  present.  Targets  are  cued  as  in  method  one  and 
reference  arrays  from  four  seeker  videos  are  loaded  sequentially  in  cor¬ 
relators  1  through  4.  This  will  require  a  maximum  of  eight  video  fields. 
The  acquisition  sensor  video  is  then  preprocessed  and  five  fields  are 
correlated  as  in  method  two  (six  fields  of  time  are  required).  Target 
assignment  and  error  signal  generation  occurs  concurrently  with  loading 
the  following  four  reference  arrays  and  therefore  requires  no  additional 

time.  After  the  correlation  surfaces  for  missile  video  five  through 

eight  are  processed  and  an  additional  four  field  delay  is  added  to  allow 
the  full  0.3  seconds  (or  18  equivalent  fields)  slewing  time,  the  first 
four  missile  LOS  are  recomputed  to  check  on  the  slewing  accuracy.  It  is 
assumed  at  this  point  that  at  least  two  of  these  four  LOS  are  within 

tracker  lock-on  accuracy  and  can  be  fired.  Error  signals  for  the  others 

will  be  generated  and  a  second  slewing  will  occur  concurrent  with  check¬ 
ing  the  slewing  accuracy  for  missiles  five  through  eight.  From  this 
point  on,  all  calculations  can  be  accomplished  within  the  missile  firing 
times  and  therefore  contribute  no  additional  time  to  the  total. 

From  the  description  of  method  three  in  the  above  paragraph,  the 
first  missile  can  be  fired  in  [3(8+6)  +  4](l/60)=  0.77  seconds  and  all 
eight  missiles  can  be  fired  in  0.77  +  7  (0.175)  =  2.00  seconds.  The  rea¬ 
son  this  method  is  faster  than  method  two  is  that  missile  LOS  slewing  of 
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four  missiles  is  occuring  simultaneously  with  calculation  of  the  LOS  of 
the  other  four  missiles  thereby  reducing  the  total  firing  time.  By  match¬ 
ing  three  missile  video  LOS  at  a  time,  six  missiles  can  be  fired  in 
[3(6+6)  +  6](l/60)  +  5  (.175)  =  1.58  seconds.  By  matching  four  missile 
video  LOS  concurrently,  four  missiles  can  be  fired  in  2 ( 8+6 ) (1/60)  +  0.3 
+  3  (.175)  =  1.29  seconds.  By  matching  only  two  missile  videos  concur¬ 
rently  then  slewing  these  two  while  the  remaining  two  are  matched,  four 
missiles  can  be  fired  in  [3(4+6)  +  8](l/60)  +  3( . 1 75 )  =  1,16  seconds. 

The  advantage  of  this  method  is  that  it  performs  faster  than  either 
method  1  or  2  and  requires  only  one-half  the  number  of  correlators  as 
method  two. 

D.  Method  Four 

Method  four  employs  three  correlators  which  operate  in  parallel. 

When  firing  on  eight  targets,  the  correlators  match  3+3+2  missile  videos. 
The  first  missile  can  be  fired  in  [(6+6)  +  (6+6)  +  (4+6 )  +  (o+6)](l/60)  = 
0.77  seconds  and  all  eight  missiles  can  be  fired  in  0.77  J  7  (.175)  = 

2.0  seconds.  Six  missiles  can  be  fired  in  [3(6+6)  +  GJ(l/60)  +  5 ( . 1 75 )  ~ 
1.58  seconds.  By  correlating  two  missile  videos  simultaneously,  four 
missiles  can  be  fired  in  1.16  seconds,  the  same  as  method  three. 

This  method  yields  identical  firing  time  results  to  method  three  but 
requires  one  less  correlator  array.  Another  advantage  of  this  method  is 
that  a  slewing  time  of  over  0.5  seconds  is  available  for  the  first  mis¬ 
sile  LOS  slewing  when  eight  targets  have  been  cued. 
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E.  Method  Five 

Method  five  employs  only  two  correlators  which  operate  in  parallel. 
The  first  missile  can  be  fired  in  5(4+6) ( 1 /60)  =  0.83  seconds  and  all 
eight  missiles  can  be  fired  1 n  0.83  +  7  (.175)  =  2.06  seconds.  The  LOS 
slewing  time  for  the  first  missile  is  0.5  seconds.  Six  missiles  can  be 
fired  in  4  (4+6 )  ( 1  /60)  +  5  ( .  1 75 )  =  1.54  seconds  and  f.>ur  missiles  can  be 
fired  in  [3(4+6)  +  8]  (1/60)  +  3  (.175)  -  1.16  seconds. 

F.  Method  Six 

In  method  six,  one  of  the  feature  matching  algorithms  discussed  in 
Chapters  2  and  3  is  used  for  multitarget  hand-off.  As  shown  in  Figure 
4-2,  targets  are  cued  in  the  same  manner  as  for  all  previous  methods. 

The  adaptive  digitizer,  prefiltering,  and  edge  detection  blocks  are  iden¬ 
tical  to  the  previous  methods.  The  segmentation  block  segments  the  scene 
into  target  and  background  and  the  feature  extractor  computes  the  feature 
vector  of  the  scene.  The  target  assignment  circuit  is  identical  to  pre¬ 
vious  methods. 

The  first  step  is  to  extract  features  from  the  missile  video  refer¬ 
ence  scenes  sequentially  and  store  in  the  feature  comparison  circuit. 
Because  the  missile  video  is  not  synchronized,  two  fields  per  reference 
are  required  for  this  step.  Features  are  then  computed  for  all  possible 
reference  size  subarrays  within  the  live  scene  (acquisition  sensor  video) 
in  real  time  and  used  by  tne  feature  comparitor  circuit  to  find  the  match 
point. 

If  one  feature  comparitor  circuit  is  used,  the  time  before  the  first 
missile  is  fired  and  the  total  time  to  fire  all  eight,  six  or  four 
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Figure  4-2.  Block  diagram  for  feature  matching  multitarget  handoff. 


missiles  are  identical  to  the  times  for  method  one.  With  eight,  four, 
three,  or  two  feature  comparitors  the  times  are  identical  to  those  for 
methods  two  through  five  respectively.  Therefore,  the  feature  matching 
methods  do  not  offer  a  time  advantage  over  the  standard  correlation 
method.  The  advantage  of  this  method  is  a  reduced  amount  of  hardware  re¬ 
quired  if  the  segmentation  and  feature  extraction  circuits  are  implement¬ 
ed  in  very  large  scale  integrated  circuits  (VL5IC)  as  the  correlator  cir¬ 
cuits  are  presently  implemented.  Less  hardware  should  be  required  be¬ 
cause  of  the  smaller  amount  of  multiplications,  additions  and  subtractions 
required  as  shown  in  previous  chapters.  If  these  VLSIC  are  not  presently 
available,  this  method  involves  more  of  a  technology  risk.  Development 
of  these  VLSIC  chips,  however,  would  be  a  step  toward  building  an  auto¬ 
matic  target,  recognizer  as  explained  in  method  seven. 

G.  Method  Seven 

Figure  4-3  is  a  block  diagram  of  the  functions  required  for  an  auto¬ 
matic  target  recognizer.  Except  for  the  target  classifier,  which  distin¬ 
guishes  a  tank  from  a  truck  or  an  APC,  etc.,  all  of  the  blocks  appear  in 
the  block  diagram  shown  in  Figure  4-2  of  a  feature  matching  multi  target 
handoff  device.  One  difference,  which  is  not  apparent  in  the  block  dia¬ 
grams,  is  that  significantly  more  feature  information  is  required  for 
target  classification  than  for  image  matching.  Also,  target  classification 
usually  requires  working  with  frames  rather  than  fields  of  viedo  in  order 
to  recognize  at  the  desired  ranges.  Developing  a  feature  matching  multi¬ 
target  handoff  device  shown  in  Figure  4-2  would  therefore  be  a  step 
toward  recognition  and  handoff  as  shown  in  Figure  4-4. 


multi  target  handoff  hardware 
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The  system  shown  in  Figure  4-4  is  initially  operated  in  the  target 
recognition  mode.  Up  to  eight  targets  are  identified  and  their  X  &  V 
positions  within  the  acquisition  sensor  FOV  are  sent  to  the  target  assign¬ 
ment  and  error  generation  circuitry  and  to  the  gunner's  monitor  where 
symbology  is  used  to  identify  targets.  Video  from  each  missile  sensor  is 
then  processed  sequentially  generating  features  of  reference  arrays  taken 
from  the  center  of  the  FOV.  It  is  anticipated  that  a  subset  of  the  fea¬ 
tures  used  for  target  recognition  will  be  used  for  image  matching.  The 
prefiltering,  edge  detection,  segmentation,  and  feature  extraction  can  all 
be  implemented  in  pipelined  very  large  scale  integrated  circuits.  There¬ 
fore,  a  small  delay  occurs  before  the  first  features  are  generated  with 
the  remainder  generated  in  real  time.  After  scaling  and  prefiltering, 
features  are  then  generated  for  all  possible  KxL  subarrays  of  the  acquisi¬ 
tion  sensor  video,  where  KxL  is  the  size  of  the  reference  arrays.  These 
features  are  then  used  to  find  the  best  match  of  each  missile  reference 
within  the  acquisition  sensor  FOV. 

With  only  one  feature  comparison  circuit,  the  time  for  firing  one, 
eight,  six,  or  four  missiles  is  id  ntical  to  the  time  for  method  one. 

With  eight,  four,  three,  or  two  feature  comparison  circuits,  the  times 
are  identical  to  methods  2,  3,  4,  or  5,  respectively.  Since  the  time 
lines  for  feature  matching  are  the  same  as  for  the  image  matching  meth¬ 
ods,  the  main  advantage  of  the  feature  matching  method  is  less  hardware 
requirement.  The  lesser  hardware  is  a  result  of  the  smaller  amount  of 
equivalent  additions  required  to  generate  the  n  correlation  surfaces. 


Tab!  e 


Method  #1 
1  Correlator 


Method  H2 
8  Correlators 


Method  #3 
4  Correlators 


Method  #4 
3  Correlators 


Method  #5 
2  Correlators 


Method  #6  &  € 
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19.  Firing 

times  for  seven  methods. 

Firing  of 

1st  Missile 
of  Eight 

8  Missiles 

6  Missiles  4  M' 

issiles 

1.2 

2.43 

1.81 

1.19 

1.07 

2.29 

1.81 

1.33 

0.77 

2.00 

1.58 

1.16 

0.77 

2.00 

1.58 

1 .16 

0.83 

2.06 

1.54 

1.16 

With  1,  8 

,  4,  3,  or  2 

feature  comparison 

circuits  times  are  equivalent  to  Methods 
1,  2,  3,  4,  or  5  respectively 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  chapter  delineates  the  major  conclusions  and  recommendations 
resulting  from  the  work  on  this  contract.  Most  of  the  conclusions  and 
recommendations  presented  in  this  chapter  have  been  given  and  justified 
within  the  first  four  chapters  of  this  report. 

A.  Conclusions 

1.  The  accuracy  and  reliability  of  several  feature  matching  digital 
image  registration  algorithms  were  investigated  through  simula¬ 
tion  using  typical  military  t/'pe  scenes.  Because  of  previous 
analysis  and  hardware  design,  a  binary  correlator  (exclusi ve-OR 
gate  implementation)  was  used  as  a  basis  for  comparison  of  the 
performance  of  feature  matchirg  algorithms  in  this  report.  When 
reference  images  were  chosen  >o  include  prominent  features  in 
the  reduced  high  resolution  FUR  images,  the  performance  of  the 
binary  correlation  algorithm  was  satisfactory.  When  each  refer¬ 
ence  image  was  chosen  from  the  center  of  the  field  of  view  of  its 
corresponding  high  resolution  image,  the  method's  performance 
was  not  as  good. 

2.  The  moments  method  and  the  method  based  on  intraset  and  interset 
distances  were  shown  not  to  be  accurate  for  the  type  imagery 
considered  in  this  work. 
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3.  The  registration  method  based  on  correlation  of  adjacent  pixels 
was  simulated  with  and  without  intensity  level  normalization  of 
images.  When  images  were  used  without  intensity  level  normaliza¬ 
tion,  the  method  failed  for  all  sixteen  pairs  of  scenes  even 
when  the  reference  images  were  chosen  to  include  the  prominent 
features  of  the  reduced  high  resolution  images.  However,  inten¬ 
sity  level  normalization  significantly  improved  the  performance 
of  the  method  indicating  the  necessity  of  some  kind  of  normali¬ 
zation  to  compensate  for  the  difference  in  the  d.c.  gains  of  the 
imaging  systems. 

4.  A  preprocessing  algorithm  which  replaces  each  pixel  by  its  value 
minus  the  average  pixel  value  of  the  KlxLl  array  centered  about 
it  compensates  for  the  difference  between  d.c.  gains  of  the  two 
imaging  systems  and  yet  allows  the  computation  of  ^  or 
V..,  .  from  V.  ..  The  performance  of  the  method  using  prepro- 

»T  *  » J  1  5 J 

cessed  images  is  comparable  to  that  of  the  binary  correlation 
algorithm. 

5.  The  performance  of  the  method  based  on  correlation  of  adjacent 
rows  and  columns  was  better  than  that  of  the  method  based  on  cor¬ 
relation  of  adjacent  pixels  in  rows  and  columns.  This  is  es¬ 
pecially  true  when  the  reference  images  have  prominent  horizontal 
and  vertical  edges.  A  schematic  for  a  possible  real-time  imple¬ 
mentation  of  the  method  of  correlation  of  adjacent  pixels  was 
given  in  Chapter  III.  This  implementation  requires  lesc  hard¬ 
ware  than  that  of  the  standard  cross  correlation  algorithm. 
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6.  In  order  to  reduce  the  hardware  requirements  further,  a  new 
method  based  on  the  summation  of  the  absolute  values  of  the  dif¬ 
ferences  between  adjacent  pixels  on  the  rows  and  on  the  columns 
was  proposed.  Simulation  results  indicate  that  the  performance 
of  this  method  is  comparable  to  the  correlation  of  adjacent 
pixels  method  discussed  in  conclusions  3,  4,  and  5. 

7.  Fast  implementations  of  the  Haar,  modified  Haar,  and  rationalized 
Haar  transforms  were  presented.  Simulations  showed  that  these 
methods  performed  comparably  with  the  correlation  of  adjacent 
pixel  methods.  The  rationalized  Haar  transform  is  easier  to  im¬ 
plement  in  hardware  and  gave  better  results  than  the  other  Haar 
methods. 

8.  Suggested  hardware  implementations  for  the  most  promising  methods 
are  given  in  Chapter  III.  The  absolute  value  of  adjacent  pixel 
difference  method  and  the  rationalized  Haar  transform  method  were 
the  simpliest  to  implement  in  hardware. 

9.  Seven  methods  for  accomplishing  automatic  bandoff  of  multiple 
targets  were  presented  in  Chapter  IV.  Table  4-1  gives  the  minimum 
time  required  to  fire  the  first  missile  and  the  minimum  time  re¬ 
quired  to  fire  eight,  six,  or  four  missiles  in  tandem  from  two 
launchers.  From  Table  4-1  it  is  clear  that  there  is  no  time  ad¬ 
vantage  to  finding  the  missile  lines  of  sight  using  one  of  the 
feature  matching  techniques.  The  advantage  of  the  feature  match¬ 
ing  tecnniques  is  less  hardware  requirements  due  to  the  smaller 
amount  of  equivalent  additions  needed  for  multiple  target  handoff. 
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The  lesser  hardware  requirement,  however,  is  conditional  on  im¬ 
plementing  the  feature  extraction  algorithms  in  LSI  curcuits. 

B.  Recommendations 

1.  If  the  desired  end  result  of  the  current  technology  program  is 
to  build  and  demonstrate  a  multiple  target  handoff  device,  it 
is  recommended  that  this  device  be  an  upgraded  version  of  the 
single  target  handoff  correlator  with  two  correlators  operating 
in  parallel.  Some  additional  circuitry  will  be  required  to  han¬ 
dle  the  sync  problems  and  some  additional  software  will  be  re¬ 
quired.  The  single  target  handoff  device  has  been  tested  and 
shown  to  work.  Upgrading  this  device  to  have  two  correlators 
operating  in  parallel  is  judged  not  to  be  a  significant  technol¬ 
ogy  risk. 

2.  If,  however,  the  objective  of  the  current  technology  program  is 
to  reduce  the  overall  helicopter  exposure  time,  some  method  of 
automatically  cuing  potential  targets  in  the  precision  pointing 
and  tracking  system  video  must  be  provided.  Moving  target  indica¬ 
tor  (MTI)  radar  or  an  automatic  target  recognizer  (or  both)  is 
needed.  It  is  well  known  that  automatic  target  recognizers  com¬ 
pare  features  of  possible  targets  with  stored  features  of  desired 
targets.  Therefore,  developing  a  multiple  target  handoff  device 
using  a  feature  matching  technique  would  be  a  step  toward  incor¬ 
porating  the  automatic  target  recognition  and  handoff  within  the 
same  hardware  device.  It  should  be  pointed  out  that  the  feature 
matching  multiple  target  handoff  technique  involves  a  greater 
technology  risk  than  the  two  parallel  correlator  technique. 
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Based  on  the  results  of  the  work  on  this  contract,  if  feature  match¬ 
ing  is  used  to  accomplish  multiple  target  handoff,  the  absolute  value  of 
adjacent  pixel  difference  method  or  the  rationalized  Haar  transform 
method  should  be  used.  To  upgrade  the  system  to  accomplish  automatic 
target  recognition  as  well  as  target  handoff  a  number  of  additional  tar¬ 
get/background  features  will  be  required. 
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